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Abstract 

The  steady  states  of  two  models  of  the  double-gyre  wind-driven  ocean  circulation  are 
studied.  The  link  between  the  steady  state  solutions  of  the  models  and  their  time-mean 
and  low-frequency  variability  is  explored  to  test  the  hypothesis  that  both  stable  and 
unstable  fixed  points  influence  shape  the  model’s  attractor  in  phase  space. 

The  steady  state  solutions  of  a  barotropic  double-gyre  ocean  model  in  which  the 
wind-stress  curl  input  of  vorticity  is  balanced  primarily  by  bottom  friction  are  studied. 
The  bifurcations  away  from  a  unique  and  stable  steady  state  are  mapped  as  a  function  of 
two  nondimensional  parameters,  {Si,5s),  which  can  be  thought  of  as  measuring  respec¬ 
tively  the  relative  importance  of  the  nonlinear  advection  and  bottom  damping  of  relative 
vorticity  to  the  advection  of  planetary  vorticity. 

A  highly  inertial  branch  characterized  by  a  circulation  with  transports  far  in  excess  of 
those  predicted  by  Sverdrup  balance  is  present  over  a  wide  range  of  parameters  including 
regions  of  parameter  space  where  other  solutions  give  more  realistic  flows.  For  the  range 
of  parameters  investigated,  in  the  limit  of  large  Reynolds  number,  Sj/Ss  oo,  the  iner¬ 
tial  branch  is  stable  and  appears  to  be  unique.  This  branch  is  anti-symmetric  with  respect 
to  the  mid-basin  latitude  like  the  prescribed  wind-stress  curl.  For  intermediate  values 
of  6i/Ss,  additional  pairs  of  mirror  image  non-symmetric  equilibria  come  into  existence. 
These  additional  equilibria  have  currents  which  redistribute  relative  vorticity  across  the 
line  of  zero  wind-stress  curl.  This  internal  redistribution  of  vorticity  prevents  the  solution 
from  developing  the  large  transports  that  are  necessary  for  the  anti-symmetric  solution 
to  achieve  a  global  vorticity  balance.  Beyond  some  critical  Reynolds  number,  the  non- 
symmetric  solutions  are  unstable  to  time-dependent  perturbations.  Time-averaged  solu¬ 
tions  in  this  parameter  regime  have  transports  comparable  in  magnitude  to  those  of  the 
non-symmetric  steady  state  branch.  Beyond  a  turning  point,  where  the  non-symmetric 
steady  state  solutions  cease  to  exist,  all  the  computed  time-dependent  model  trajectories 
converge  to  the  anti-symmetric  inertial  runaway  solution.  The  internal  compensation 
mechanism  which  acts  through  explicitly  simulated  eddies  is  itself  dependent  explicit 
dissipation  parameter. 
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Using  the  reduced-gravity  quasigeostrophic  model  an  investigation  of  the  link  between 
the  steady  state  solutions  and  the  model’s  low-frequency  variability  is  conducted.  If  the 
wind-stress  curl  is  kept  anti-symmetric,  successive  pairs  of  non-symmetric  equilibria  come 
into  existence  via  symmetry-breaking  pitchfork  bifurcations  as  the  model’s  biharmonic 
viscosity  is  reduced.  Succesive  pairs  of  mirror  image  equilibria  have  an  additional  half 
meander  in  the  jet.  The  distinct  energy  levels  of  the  steady  state  solutions  can  be  under¬ 
stood  in  part  by  there  different  inter-gyre  fluxes  of  vorticity.  Those  solutions  with  weak 
inter-gyre  fluxes  of  vorticity  have  large  and  energetic  recirculation  cells  which  remove 
excess  vorticity  through  bottom  friction.  Those  solutions  with  strong  inter-gyre  fluxes  of 
vorticity  have  much  smaller  and  less  energetic  recirculation  cells. 

A  significant  fraction  of  the  variance  (30%)  of  the  interface  height  anomaly  can  be 
accounted  by  four  coherent  structures  which  point  away  from  the  time-mean  state  and 
towards  four  steady  state  solutions  in  phase  space.  After  removing  the  variance  which 
projects  onto  the  four  modes,  the  remaining  variance  is  reduced  predominantly  at  low- 
frequencies,  showing  that  these  modes  are  linked  to  the  low-frequency  variability  of  the 
model.  Furthermore,  the  time-averaged  flow  fields  within  distinct  energy  ranges  show 
distinct  patterns  which  are  in  turn  similar  to  the  distinct  steady  state  solutions. 
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2.1  Parameter  chart  for  5h  =  0.04,  showing  the  bifurcations  of  the  anti¬ 

symmetric  equilibrium  (A).  Dashed  lines  indicate  pitchfork  bifurcations 
and  dotted  lines  indicate  Hopf  bifurcations.  The  anti-symmetric  equi¬ 
librium  exist  for  the  full  range  of  parameters  shown  in  the  chart.  It  is 
unstable  in  the  shaded  region,  and  stable  in  the  unshaded  region.  The 
lower  dashed  lobe,  indicates  where  non-symmetric  equilibria  of  type  N1 
(lower  part  of  lower  dashed  lobe),  and  of  type  N2  (upper  part  of  lower 
dashed  lobe)  bifurcate  from  the  anti-symmetric  equilibrium  via  a  symme¬ 
try  breaking  pitchfork  bifurcation.  The  upper  dashed  lobe  indicates  where 
non-symmetric  equilibria  of  type  N3  bifurcate  from  the  anti-symmetric 
branch  via  a  pitchfork  bifurcation .  36 

2.2  Contour  plots  of  xf)  (top  row)  and  q  (bottom  row)  for  the  anti-symmetric 

equilibrium  (type  A)  with  5s  =  0.01,  5h  =  0.04,  and  5j  as  indicated.  The 
dashed  lines  indicate  the  negative  contours  and  the  solid  lines  indicate  the 
positive  contours.  The  thick  solid  line  indicates  the  zero  contour.  The 
contour  interval  is  0.03  for  xj;  and  0.02  for  9 .  37 
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2.3  Parameter  chart  for  =  0.04,  showing  the  bifurcations  of  the  non- 

symmetric  equilibria  of  type  Nl.  Solid  lines  denote  saddle-node  bifurca¬ 
tions,  dashed  lines  indicate  pitchfork  bifurcations,  and  dotted  lines  indicate 
Hopf  bifurcations.  Equilibria  of  type  Nl  can  be  traced  continuously  for 
increasing  values  of  Sj  from  the  dashed  curve  (where  they  bifurcate  from 
the  anti-symmetric  equilibrium),  up  to  the  saddle-node  curve  labeled  SN 
or  SNl  depending  on  whether  8s  is  to  the  right  or  left  of  the  cusp  point 
labeled  CPI .  40 

2.4  Typical  ip  and  q  fields  for  the  branches  A,N1  and  N2.  Sj  =  0.1001,  Ss  = 

0.02  and  Sh  =  0.04.  The  contour  interval  is  0.3  for  ip  and  0.2  for  q.  The 
negative  contours  are  dashed,  and  the  zero  contour  is  the  thick  one.  ...  41 

2.5  Bifurcation  plot  of  ipmax  +  -vs-  5/,  indicating  the  emergence  of  non- 
symmetric  equilibria  Nl  and  N2  and  N3  via  pitchfork  bifurcations  at 
5 1  =  5pi,  5p2,  and  ^P3.  A  fourth  pitchfork  bifurcation  point  at  8j  =  Sp4 
marks  the  disappearance  of  the  non-symmetric  equilibria  of  type  N3.  Anti¬ 
symmetric  equilibria  lie  on  the  horizontal  line  in  the  center  of  the  figure. 
Saddle-node  bifurcation  points  at  5/  -  ^51,  ^52,  and  833,  mark  the  merging 
of  equilibria  of  type  Nl  with  Nli,  Nli  with  NI2,  and  of  type  NI2  with  N2. 

The  solid  curves  indicate  equilibria  which  are  stable  and  the  dashed  lines 
indicate  equilibria  which  are  unstable.  {83  =  0.01,  8h  =  0.04) .  42 

2.6  Contour  plots  of  ip  for  81  =  0.0316,  83  =  0.001  and  8h  =  0.0313.  Steady 

states  of  type  A,  Nl,  and  N2  are  presented  along  with  the  time- averages 
stream-function  field  (PBAR).  The  solid  lines  denote  positive  contours 
and  the  dashed  lines  denote  negative  contours.  Contour  intervals  are  as 
indicated .  43 
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2.7  Parameter  chart  for  5ji  —  0.04,  showing  the  bifurcations  of  the  non- 
symmetric  equilibrium  (N2).  Solid  lines  denote  saddle-node  bifurcations 
and  dashed  lines  indicate  pitchfork  bifurcations.  Equilibria  of  type  N2 
exist  only  in  the  shaded  region  of  parameter  space  where  it  is  also  al¬ 
ways  unstable.  Non-symmetric  equilibria  of  type  N2  come  into  existence 
via  a  pitchfork  bifurcation  of  the  anti-symmetric  flow,  and  can  be  traced 
continuously  as  5/  is  increased  up  to  the  solid  curve  labeled  SN  where  it 
experiences  a  saddle-node  bifurcation,  merges  with  equilibria  of  type  Nl, 

and  ceases  to  exist .  45 

2.8  Typical  and  q  fields  for  the  branches  Al  and  N3.  5j  =  0.1127,  5s  =  0.01 
and  5h  =  0.04.  The  contour  interval  is  0.4  for  if)  and  0.3  for  q.  The 
negative  contours  are  dashed,  and  the  zero  contour  is  the  thick  one.  The 
non-symmetric  branch  N3  has  its  inter-gyre  boundary  shifted  northward 

and  tilted  from  west  to  east .  47 

2.9  Typical  V’  and  q  fields  for  the  branches  Nl,  Nli,  and  NI2,  for  5i  =  0.1, 

5s  =  0.02  and  5h  =  0.04.  The  contour  interval  is  0.2  for  ^  and  0.2  for 
q.  The  negative  contours  are  dashed,  and  the  zero  contour  is  the  thick 
one.  The  difference  between  the  three  equilibrium  branches  is  essentially 
restricted  to  the  spatial  arrangement  of  the  multiple  closed  circulation  cells.  48 

2.10  Contour  plots  of  ^p  (top  row)  and  q  (bottom  row)  for  the  anti-symmetric 
equilibrium  (type  A)  with  5s  =  0.04  and  5h  =  0.01  and  5/  as  indicated. 

The  dashed  lines  indicate  the  negative  contours  and  the  solid  lines  indicate 
the  positive  contours.  The  thick  solid  line  indicates  the  zero  contour.  C.I. 

=  0.02  for  Ip  and  C.I.  =  0.02  for  o' .  50 
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2.11  Contour  plots  of  ^1)  and  q  for  the  single-gyre  case.  The  no-flux  of  vorticity 
across  the  northern  wall  allows  the  northern  jet  to  penetrate  straight  across 
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2.12  Bifurcation  plot  of  maxip  +  min^  -vs-  Sj,  indicating  the  emergence  of  non- 
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2.13  Contour  plots  of  ij)  (top  row)  and  q  (bottom  row)  around  the  nose  point 
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2.14  Parameter  chart  for  5h  =  0.04  showing  an  overlay  of  the  bifurcations 
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Chapter  1 


Introduction 


Numerical  ocean  models  are  an  indispensable  tool  for  understanding  the  climate  system 
and  possibly  for  predicting  climate  change.  Ocean  models  are  not  only  used  in  conjunc¬ 
tion  with  observations  to  estimate  the  current  state  of  the  oceans,  but  also  to  estimate 
the  state  of  the  ocean  under  different  mechanical  and  thermodynamical  forcing.  These 
models  depend  on  boundary  conditions  and  sub-grid-scale  parameterizations  that  are 
poorly  known  from  observations.  For  climate  studies,  the  time  evolution  of  ocean  models 
over  hundreds  to  thousands  of  years  is  of  paramount  importance.  This  makes  the  choice 
of  suitable  parameterizations  of  dissipation  rather  crucial,  since  dissipative  forces,  no 
matter  how  small,  have  enough  time  to  become  important. 

For  the  wind-driven  ocean  circulation,  the  sub-grid-scale  parameterization  of  mix¬ 
ing  processes  provides  an  explicit  dissipation  term  in  the  governing  equation.  Pedlosky 
(1996),  reviews  the  role  played  by  dissipation  in  theories  of  the  wind-driven  circulation 
within  the  context  of  homogeneous  models.  In  all  cases  that  have  been  studied,  the  ex¬ 
plicit  frictional  dissipation  is  responsible  for  balancing  the  continuous  input  of  vorticity 
by  the  action  of  the  wind-stress.  The  hope  that  as  the  boundary-layer  Reynolds  number 
is  increased,  the  total  circulation  would  become  independent  of  the  particular  frictional 
model  adopted  has  been  disappointed.  As  demonstrated  by  lerley  and  Sheremet  (1995), 
in  a  single-gyre  model  with  free-slip  boundary  conditions,  a  steady  basin-filling  inertial 
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gyre  with  velocities  far  in  excess  of  those  predicted  by  Sverdrup  balance  is  the  only 
solution  for  sufficiently  large  boundary-layer  Reynolds  number. 

Of  particular  interest  to  our  study  is  the  double-gyre  simulation  by  Marshall  (1984). 
In  this  model,  a  source  of  negative  vorticity  input  by  the  curl  of  the  wind-stress  in 
the  southern  region  of  the  basin  is  balanced  by  a  positive  source  of  vorticity  of  equal 
magnitude  in  the  northern  region.  Marshall’s,  study  stands  out  because  it  appears  to 
provide  an  example  in  which  time-dependent  eddies  prevent  the  development  of  the 
inertial  runaway  solution.  In  this  simulation,  the  amount  of  negative  vorticity  put  in 
by  the  wind  in  the  sub-tropical  gyre  is  roughly  balanced  by  the  eddy  flux  of  negative 
vorticity  from  the  southern  gyre  to  the  northern  gyre,  thereby  eliminating  the  need 
for  the  vorticity  input  by  the  wind  to  be  eliminated  in  the  western  boundary-layer.  It 
is  important  to  point  out,  however,  as  Pedlosky  (1996),  emphasizes,  that  this  internal 
compensation  mechanism  can  only  apply  for  the  singular  case  in  which  there  is  no  net 
input  of  vorticity  by  the  curl  of  the  wind-stress  over  the  entire  domain  —  any  imbalance 
must  be  removed  by  the  explicit  dissipation. 

Despite  its  limited  applicability  to  the  real  ocean  case,  the  internal  compensation 
mechanism  deserves  further  study  since  its  action  appears  to  make  the  time-averaged  so¬ 
lutions  independent  of  the  explicit  dissipation  parameterization.  In  particular,  one  would 
want  to  know  which  type  of  instabilities  allow  the  internal  compensation  mechanism  to 
act,  and  whether  inertial  runaway  can  be  truly  avoided  as  the  boundary-layer  Reynolds 
number  is  increased.  The  first  motivation  for  this  study  is  to  address  these  issues. 

To  this  end,  we  investigate  the  steady  state  solutions  and  their  stability  for  a  large 
range  of  parameters.  The  techniques  of  numerical  bifurcation  theory  (Seydel  1994)  are 
used  to  unravel  the  bifurcation  structure  of  the  steady  state  equilibria  of  the  ocean  model 
introduced  by  Marshall  (1984).  The  model’s  interesting  novelty  is  that  the  boundary  con¬ 
ditions  are  such  that  no  relative  vorticity  flux  is  allowed  through  the  basin  walls,  despite 
the  fact  that  the  model  has  both  bottom  friction  and  lateral  diffusion.  The  boundary  con¬ 
ditions  are  therefore  dynamically  equivalent  to  those  of  a  model  having  bottom  friction 
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alone.  The  model  nonetheless  retains  a  lateral  diffusion  term  which  prevents  the  devel¬ 
opment  of  discontinuities  in  the  relative  vorticity  field.  This  formulation  of  the  boundary 
condition  for  the  eddy  diffusion  term  is  based  on  the  observation  that  geostrophic  eddies 
act  only  to  redistribute  vorticity  laterally  and  not  as  a  sink  of  vorticity  through  the  basin 
walls.  The  model  is  essentially  the  same  as  that  used  in  a  study  of  steady  state  solutions 
by  Cessi  and  lerley  (1995),  but  differs  in  the  choice  of  dissipation  operator  and  boundary 
conditions;  Cessi  and  lerley  (1995),  used  Munk  type  lateral  diffusion  with  free-slip  at  the 
eastern  and  western  walls  and  periodic  boundary  conditions  at  the  northern  and  south¬ 
ern  walls.  A  discussion  of  this  super-slip  boundary  condition,  as  well  as  other  choices  of 
boundary  conditions  can  be  found  in  Pedlosky  (1996). 

By  studying  the  steady  state  solutions,  we  will  discover  in  which  region  of  parameter 
space  the  circulation  retains  the  Sverdrup  balance  as  part  of  the  solution  and  in  which 
part  of  parameter  space  the  circulation  is  of  the  inertial  runaway  type.  By  studying  the 
bifurcation  structure  of  the  steady  state  solutions,  we  will  map  out  where  qualitative 
changes  in  the  nature  of  the  solution  occur,  and  thus  carve  out  regions  of  parameter 
space  where  the  internal  compensation  mechanism  can  act. 

The  second  motivation  for  this  study  deals  with  the  issues  of  low-frequency  variability 
and  multiple  equilibria.  As  pointed  out  by  Jiang  et  al.  (1995),  the  oceans’  western 
boundary  currents  offer  clear  examples  of  low-frequency  variability  in  the  wind-driven 
ocean  circulation.  Some  examples  are  provided  by  the  path  of  the  Kuroshio  alternating 
between  a  large  and  a  small  meander  state  with  a  period  of  several  years  (Taft  1972),  by 
the  latitude  of  separation  of  the  Brazil/Malvinas  current  system  varying  on  inter-annual 
time-scales,  (Olson  et  al.  1988),  and  by  the  mean  position  of  the  Gulf  stream  that  varies 
inter-annually  (Brown  and  Evans  1987).  Inter-annual  and  longer  time-scale  variability  is 
a  possible  manifestation  of  the  nonlinearity  of  the  wind-driven  circulation  since  it  cannot 
be  accounted  for  by  the  seasonal  cycle  of  the  forcing.  Other  mechanisms  are  however 
possible.  Frankignoul  and  Hasselman  (1977),  have  shown  that  because  of  the  relatively 
slow  response  time  of  the  ocean  compared  to  the  atmosphere,  the  ocean  can  integrate  the 
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high  frequency  atmospheric  forcing  leading  to  enhanced  variability  at  longer  time-scales. 
Another  possibility  is  that  some  of  the  low  frequency  variability  in  the  climate  system 
can  be  understood  in  terms  of  coupled  modes  between  the  ocean  and  atmosphere.  In  this 
scenario,  temperature  anomalies  advected  around  the  oceanic  gyres  with  a  time-scale  of 
a  decade  or  so,  couple  to  the  atmosphere  to  produce  wind  fields  which  in  turn  re-enforce 
the  anomalies.  Finally,  another  possibility  is  that  the  variability  is  due  to  the  internal 
dynamics  of  ocean  currents.  As  we  will  discuss  below,  adiabatic  ocean  models  forced  by 
steady  winds  are  quite  capable  of  producing  the  types  of  variability  offered  by  the  above 
examples. 

Because  of  the  nonlinearity  of  the  governing  equations,  ocean  models  can  exhibit  a  rich 
variety  of  dynamical  behaviors,  including  multiple  equilibria,  self  sustained  oscillations 
and  chaos.  The  ideas  of  dynamical  systems  theory  can  be  used  to  investigate  the  intrinsic 
variability  of  these  models.  Dynamical  systems  theory  has  been  applied  extensively  to 
the  thermohaline  circulation  because  of  its  clear  role  in  the  earth’s  climate  system  — 
it  transports  large  amounts  of  heat  poleward.  The  gyre  dynamics  associated  with  the 
wind-driven  circulation  has  received  much  less  attention  from  the  point  of  view  of  low- 
frequency  variability  and  multiple  equilibria.  However,  t  is  surely  not  less  important. 
For  example,  the  mid-latitude  gyres  of  the  North  Pacific  are  probably  the  main  agents 
of  poleward  heat  transport  in  that  ocean,  and  in  the  North  Atlantic  the  modeling  study 
of  Spall  (1996a),  (1996b),  has  demonstrated  how  the  surface  wind-driven  circulation  can 
be  coupled  to  the  deep  western  boundary  current  and  thus  affect  the  strength  of  the 
thermohaline  circulation. 

Recently,  several  studies  have  introduced  the  concepts  of  attractors  and  fixed  points 
to  help  characterize  the  behavior  of  wind-driven  ocean  models.  For  dissipative  dynamical 
systems,  all  solutions  converge,  as  t  oo,  to  a  complicated  set  called  the  global  attractor, 
which  may  be  fractal.  This  set  is,  in  general,  finite  dimensional,  and  embodies  the  long 
time  evolution  of  the  model,  including  turbulent  states.  Although  it  is  not  yet  numerically 
feasible  to  fully  map  out  the  global  attractor  of  general  circulation  ocean  models,  it 
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is  possible  to  find  their  fixed  points  .  These  fixed  points  lie  in  the  global  attractor 
and  can  be  useful  in  characterizing  the  long  time  evolution  of  the  model.  Speich  et 
al.  (1995)  extended  the  work  of  Jiang  (1995)  and  mapped  the  successive  bifurcations  of 
a  wind-driven  shallow-water  model.  They  found  multiple  equilibria  for  realistic  values 
of  the  forcing  and  dissipation  parameters.  These  equilibria  are  the  result  of  nonlinear 
interactions  between  the  two  main  recirculation  cells  that  flank  the  sides  of  the  inter¬ 
gyre  jet.  The  equilibria  come  into  existence  via  an  imperfect  or  perturbed  pitchfork 
bifurcation  of  the  nearly  anti-symmetric  solution  ^  (see  Appendix  A).  They  also  made 
some  comparison  with  observations  of  the  Gulf  Stream  and  Kuroshio  Extensions  and 
concluded  that  the  internal  variability  of  their  simulated  ocean  could  be  an  important 
factor  in  the  observed  inter-annual  variability  of  these  currents. 

Dijkstra  and  Katsman  (1997)  explored  the  successive  transitions  from  a  unique  steady 
state  solution  to  time-dependence  of  the  doulbe-gyre  circulation  in  a  reduced  gravity  and 
in  a  2  layer  quasigeostrophic  model  as  the  lateral  friction  was  decreased.  For  the  reduced 
gravity  model  they  identified  two  classes  of  modes  which  lead  to  oscillatory  instabilities. 
The  first  class  with  a  monthly  to  inter-monthly  period  originates  from  the  ocean  basin 
modes  of  the  inviscid  theory  (Pedlosky  1987).  Another  class  of  modes  which  they  call 
gyre-modes  exists  only  because  of  the  presence  of  the  gyres.  These  modes  have  an  inter¬ 
annual  period  of  the  order  of  the  advective  time  scale  of  the  gyre  circulation.  Both  modes 
were  found  to  become  unstable  for  values  of  lateral  friction  lower  than  that  at  which  the 
symmetry  breaking  pitchfork  bifurcation  occurred.  This  should  be  contrasted  with  the 
2-layer  model  where  the  linear  stability  boundary  of  the  anti-symmetric  solution  consists 
of  a  Hopf  bifurcation,  and  not  the  pitchfork  bifurcation  (Dijkstra  and  Katsman  1997). 
The  Hopf  bifurcation  is  the  result  of  baroclinic  instability  and  the  excited  modes  have  a 
period  with  the  annual  time-scale. 

Meacham  (submitted  1997)  studied  the  origin  of  low-frequency  variability  in  a  ho- 

^For  these  simulations  the  wind-stress  curl  profile  is  anti-symmetric,  but  the  shallow  water  governing 
equations  do  not  have  the  meridional  symmetry  property  of  the  QG  equations. 
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mogeneous  ocean  model  forced  by  steady  wind-stress.  He  found  that  irregular  large 
amplitude  vacillations  associated  with  the  inertial  recirculation  gyres  characterize  the 
behavior  of  the  model  when  the  forcing  is  sufficiently  strong.  In  addition,  he  found  three 
limit-cycles  which  could  be  followed  over  ranges  of  parameters  where  they  were  stable 
and  conjectured  that  these  simple  solutions  continued  to  exist,  though  unstable,  for  lower 
values  of  the  dissipation  parameter.  These  limit  cycles  come  into  existence  via  a  Hopf 
bifurcation  for  each  of  the  three  steady  state  solutions,  one  anti-symmetric  and  the  other 
two  from  a  pair  of  non-symmetric  equilibria.  He  also  showed  that  the  large  amplitude  ir¬ 
regular  vacillation  is  most  likely  a  cycling  between  the  neighborhood  of  the  non-symmetric 
limit  cycles  and  the  low  energy  anti-symmetric  limit  cycle.  The  type  of  variability  that 
results  from  such  oscillations  depends  on  the  behavior  of  phase  space  along  the  whole 
trajectory  and  is  the  result  of  what  is  known  as  a  global  bifurcation.  In  general  the  fixed 
point  solutions  have  associated  with  them  a  stable  and  an  unstable  manifold.  Almost 
all  trajectories  which  begin  in  a  neighborhood  of  the  fixed  point  are  eventually  expelled 
along  the  unstable  manifold.  These  trajectories  can,  however,  sometimes  return  to  the 
neighborhood  of  the  fixed  point  by  following  a  trajectory  close  to  the  stable  manifold.  In 
this  sense  the  unstable  fixed  points  act  to  “steer”  the  trajectory  of  the  time  dependent 
models  (Legras  and  Ghil  1985),  thus  providing  a  mechanism  for  recurrent  and  persistent 
dynamical  regimes.  Global  bifurcations  happen  when  the  unstable  manifold  intersects 
the  stable  manifold  of  a  fixed  point  thus  giving  rise  to  a  large  amplitude  oscillation.  This 
should  be  contrasted  with  the  oscillations  associated  with  Hopf  bifurcation  points  which 
give  rise  to  oscillations  whose  amplitude  grow  as  where  Uf.  is  the  critical  value  of 

the  parameter  at  the  bifurcation  point.  Hopf  bifurcations  are  local  in  that  the  resulting 
oscillation  depends  only  on  the  local  behavior  of  phase  space  and  can  be  understood  from 
the  model  linearized  about  the  bifurcation  point.  Global  bifurcations  on  the  other  hand 
have  oscillations  which  are  of  finite  amplitude,  and  can  not  be  understood  only  in  terms 
of  the  linear  behavior  of  the  model  near  the  fixed  point. 

The  idea  that  some  simple  though  unstable  solutions  can  strongly  influence  the  shape 
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of  the  system  attractor  is  not  a  new  one,  and  has  been  explored  to  some  extent  in 
atmospheric  models,  in  the  context  of  preferred  weather  regimes.  It  has  been  shown  in 
simple  atmospheric  models  that  the  time-dependent  trajectory  can  spend  a  considerable 
amount  of  time  in  the  neighborhood  of  simple  steady  state  solutions  even  if  they  are 
unstable  (Legras  and  Ghil  1985;  Branstator  and  Opsteegh  1989). 

Preferred  dynamical  regimes  are  also  present  in  ocean  models.  Using  a  primitive 
equation  model  of  the  Gulf  Stream/Deep  Western  Boundary  Current  crossover.  Spall 
(1996b),  found  low-frequency  variability  associated  with  the  transition  between  two  pre¬ 
ferred  dynamical  states,  which  included  a  high  energy  state  with  the  Gulf  stream  exten¬ 
sion  penetrating  deep  into  the  basin  and  a  low  energy  state  with  a  weakly  penetrating 
Gulf  Stream  extension.  The  mechanism  controlling  the  transition  between  the  two  states 
involved  the  interaction  between  the  surface  wind-driven  currents  and  the  deep  western 
boundary  current. 

In  a  simpler  reduced  gravity  QG  model  of  the  double-gyre  circulation,  McCalpin  and 
Haidvogel  (1996),  also  found  low-frequency  variability.  It  was  associated  with  irregular 
transitions  among  several  preferred  dynamical  regimes,  including  a  high  energy  state  with 
a  jet  penetrating  deep  into  the  basin,  a  low  energy  state  with  a  weakly  penetrating  jet, 
and  an  intermediate  energy  state  with  intermediate  jet  penetration.  Despite  the  differ¬ 
ences  between  the  models  of  Spall  and  McCalpin  and  Haidvogel,  the  preferred  dynamical 
regimes  were  similar  in  both  studies.  The  jet  penetration  scale  and  the  intensity  of  the 
eddy  energy  field  varied  among  states  in  a  similar  fashion  for  both  models.  Even  though 
the  mechanism  for  the  transition  between  states  is  different,  the  similar  nature  of  the 
preferred  regimes  suggest  that  the  existence  of  these  regimes  might  be  the  result  of  the 
dynamics  of  the  wind-driven  circulation  alone.  It  further  suggests  the  possibility  that  the 
existence  of  multiple  regimes  might  persist  through  a  hierarchy  of  models,  from  simple 
QG  dynamics  to  progressively  more  sophisticated  dynamics  and  perhaps  to  fully  coupled 
ocean  atmosphere  climate  models. 

The  above  considerations  make  it  crucial  that  the  long  time  evolution  (t  oo)  of 
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ocean  models  of  varying  complexity  and  with  differing  sub-grid-scale  parameterizations 
and  boundary  conditions  be  investigated  in  a  systematic  way  for  a  wide  range  of  parame¬ 
ter  values.  As  a  contribution  towards  this  goal,  this  thesis  presents  a  study  of  the  steady 
state  solutions  of  homogeneous  models  of  the  non-linear  wind-driven  ocean  circulation. 
Chapter  2  consists  of  a  study  of  the  the  fixed  points  and  of  the  fixed  points’  stability 
properties  for  a  barotropic  quasigeostrophic  model  with  super-slip  boundary  condition. 
Chapter  3  focuses  on  intrinsic  low-frequency  variability  within  the  context  of  a  reduced 
gravity  quasigeostrophic  model  with  free-slip  boundary  conditions.  The  steady  state  solu¬ 
tions  are  independent  of  the  choice  of  deformation  radius  since  the  radius  of  deformation 
enters  the  governing  equation  only  in  a  term  with  a  time  derivative.  Because  of  this,  the 
fixed  points  of  the  barotropic  model  are  also  fixed  points  of  the  reduced  gravity  model, 
and  vice  versa.  The  stability  of  the  fixed  points  does  however  depend  on  the  choice  of 
deformation  radius.  In  Chapter  4,  we  review  the  thesis’  major  results  and  discuss  possible 
future  research  directions.  Appendix  A  reviews  the  some  terminology  from  bifurcation 
theory. 
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Chapter  2 


Barotropic  Model 

2.1  Homogeneous  Models  of  the  Wind-Driven  Cir¬ 
culation:  Review 

It  is  useful  to  review  the  theory  of  the  wind-driven  ocean  circulation  within  the  context 
of  homogeneous  models  before  proceeding.  A  more  complete  discussion  can  be  found 
in  Pedlosky  (1996).  Scaling  analysis  suggests  that  over  most  of  the  mid-latitude  ocean 
basins,  inertial  and  frictional  terms  in  the  vorticity  equation  can  be  neglected  in  favor  of  a 
balance  between  the  advection  of  planetary  vorticty  gradients  and  the  input  of  vorticity  by 
the  curl  of  the  wind-stress.  This  linear  vorticity  balance,  also  called  the  Sverdrup  relation 
cannot,  however,  hold  over  the  whole  basin,  since  dissipative  effects  must  be  important  in 
some  part  of  the  basin  if  the  flow  is  to  remain  bounded  in  time.  Through  the  inclusion  of  a 
boundary  layer  along  the  western  edge  of  the  basin,  disspative  effects  can  be  incorporated 
into  the  solution  without  the  need  of  abandoning  the  Sverdrup  solution  over  the  rest  of 
the  basin.  But  when  the  resulting  solution  is  substituted  into  the  discarded  inertial 
terms,  it  is  discovered  that  in  the  boundary  layer,  these  terms  are  as  important  as  others 
in  the  vorticity  equation.  Therefore,  the  inertial  terms  must  be  retained  in  the  western 
boundary  layer  if  nowhere  else.  It  turns  out  that  where  the  flow  field  of  the  Sverdrup 


30 


solution  is  westward,  the  effect  of  the  inertial  terms  can  be  added  successfully  to  the 
solution  without  destroying  the  boundary  layer  nature  of  the  solution.  Where  the  flow 
field  is  eastward,  however,  the  solution  loses  its  boundary  layer  character.  In  numerical 
simulations,  these  regions  of  eastward  flow  are  the  site  where  recirculation  gyres  form. 
Furthermore,  as  the  boundary  layer  Reynolds  number  (measuring  the  relative  strength 
of  inertial  and  frictional  terms)  increasese,  these  recirculation  gyres  extend  outwards  to 
fill  the  entire  basin  and  completely  destroy  the  Sverdrup  solution. 


2.2  Model  Formulation 


As  mentioned  in  the  introduction,  the  model  configuration  is  the  same  as  that  used  in 
Marshall  (1984).  The  governing  equation,  in  nondimensional  form,  is  the  barotropic 
vorticity  equation  with  bottom  friction  and  biharmonic  lateral  diffusion, 

^  +  XT-SsC-SfiV^inV,  (2.1) 

where 

C  =  V^ip  and  q  =  S](  +  y,  (2.2) 


are  the  relative  and  potential  vorticities. 

The  dimensionless  parameters  in  the  problem  are: 
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the  inertial,  Stommel  and  diffusive  layer  thicknesses  scaled  by  the  width  of  the  basin 
L.  The  scales  which  lead  to  the  nondimensional  parameters  and  that  must  be  used  to 
reconstruct  the  dimensional  variables  are  the  following: 
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The  domain  of  integration  is  a  rectangular  basin  given  by 


(2.4) 


T>  =  {(a:,  y)  |  0  <  a;  <  1  and  —  1  <  y  <  1}. 


(2.5) 
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The  boundary  conditions  are  the  following 


^  =  0,  on  dT>,  (2.6) 

VC  •  n  =  0,  on  dT>,  (2.7) 

V(V^C)  •  n  =  0  on  (2.8) 

Note  that  the  integral  of  the  lateral  diffusive  term  over  the  entire  basin  vanishes  because 
of  the  no-flux  boundary  conditions  in  Eq.  7,  so  that  no  net  source  or  sink  of  vorticity 
is  introduced  by  the  lateral  diffusion.  The  source  of  vorticity  due  to  the  curl  of  the 
wind-stress  is  given  by 

V  X  r  =  7rsin(7ry).  (2.9) 

As  Veronis  (1966)  discusses,  the  fact  that  many  complicating  physical  processes  are 
assumed  out  of  the  system  means  that  we  can  not  think  of  such  a  simple  model  as  the 
first  term  in  a  sequence  that  converges  to  the  “real”  ocean.  Rather,  the  utility  of  such 
a  model  is  that  it  can  be  used  to  check  and  build  our  intuition  about  the  behavior  of 
oceanic  models.  In  this  respect,  it  is  better  thought  of  as  being  at  the  base  of  a  hierarchy 
of  models  which  have  successive  levels  of  sophistication  and  realism.  Investigating  the 
behavior  of  the  model  as  parameters  are  allowed  to  tend  to  various  limits  is  fundamental 
in  characterizing  its  behavior.  Thus,  the  goal  of  this  study  is  to  consider  a  wide  range  of 
parameter  (J/,  6s,  6h)  values. 

It  is  useful  to  consider  some  typical  values  for  the  model’s  parameters  which  give 
results  reproducing  the  major  gross  features  of  the  wind-driven  circulation.  Following 
Marshall  (1984),  the  inertial  layer  thickness,  6i  =  \/l  x  10“^  =  0.0316  corresponds  to  an 
ocean  in  which  r  =  10“^  N  m"^,  p  =  10^  kg  m“^,  /?  =  2  x  10"^^  m-^s“\  =  5  x  10^  m, 

L  =  10®  m;  the  nondimensional  Stommel  layer  thickness  6s  =  10“^  corresponds  to  /  = 
10“^  s~\  A„  =  1  X  lO""*  m^  s“\  and  the  horizontal  hyper-diffusive  thickness  6h  =  0.0313, 
corresponds  to  Ah  =  6  x  10^^  x  m^  s“^.  The  time  scale  is  given  by  T  =  l/PL,  so  that  one 
can  relate  6s  to  a  damping  time-scale  oi  a  —  l/{6sPL)  =  58  days.  The  lateral  diffusion 
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parameter  is  usually  chosen  to  be  the  smallest  possible  value  that  will  keep  the  vorticity 
fields  free  of  grid  scale  oscillations,  but  is  often  justified  as  an  effective  eddy  viscosity. 

Many  physical  processes  contribute  to  the  effective  dissipation,  making  it  difficult  to 
estimate  from  observation.  Most  of  these  processes  cannot  be  directly  simulated  in  a 
barotropic  QG  model  and  are  lumped  in  the  effective  eddy  diffusivity,  while  others  like 
the  internal  compensation  mechanism,  can  be  simulated  provided  the  model  resolution  is 
sufficient  to  allow  eddies.  The  model  behavior  is  relevant  even  where  it  does  not  produce 
realistic  flows,  since  it  helps  to  understand  how  processes  that  are  being  simulated  depend 
on  those  processes  parameterized  by  eddy  diffusivity.  The  method  of  solution  is  discussed 
in  detail  in  Appendix  B. 

2.3  Multiple  Equilibria 

2.3.1  Bifurcation  Structure  as  a  Function  of  6i  and  6s 

Using  a  continuation  algorithm  for  finding  both  steady  state  solutions  and  the  corre¬ 
sponding  least  stable  eigenmode  (or  one  of  the  unstable  eigenmodes  if  the  solution  is 
unstable),  the  bifurcations  of  the  steady  state  equilibria  were  mapped  as  a  function  of 
the  nondimensional  parameters  Sj  and  5s.  All  of  the  solutions  presented  in  this  section 
used  Nx  =  33  and  Ny  =  65  grid  points  in  the  x  and  y  directions  respectively,  with 
uniform  grid  spacing  of  dx  =  dy  =  0.03.  The  lateral  diffusivity  was  fixed  at  =  0.04 
which  is  slightly  larger  than  the  grid  spacing.  Section  2.3.7  discusses  the  sensitivity  of 
the  solutions  to  the  magnitude  of  5//. 

For  parameter  values  in  the  range  0.01  <  ^5  <  0.1  and  0  <  <5/  <  0.4,  up  to  6  different 
types  of  equilibria  were  found.  Each  will  be  discussed  in  turn  in  the  next  subsections. 
Before  proceeding  it  is  useful  to  give  a  review  of  the  multiple  equilibria  results  of  Cessi 
and  lerley  (1995),  followed  by  a  brief  overview  of  the  multiple  equilibria  found  in  this 
study. 

It  is  important  to  recall  that  the  model  formulation  used  by  Cessi  and  lerley,  as  well  as 
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the  one  used  here,  satisfies  the  symmetry  property  ^  >  —ip,  y  — >  —y.  Equilibrium  states 
that  satisfy  the  above  symmetry  property  are  said  to  be  anti-symmetric,  and  those  that 
do  not  are  said  to  be  non-symmetric.  Because  of  this  symmetry  property,  non-symmetric 
equilibria  always  come  in  pairs  which  are  related  to  each  other  by  ^  y  ->  —y. 

Cessi  and  lerley  (1995),  identified  5  different  types  of  equilibria  in  a  parameter  space 
defined  by  (<5/,  5a/),  where  Sm  is  the  nondimensional  Munk  boundary-layer  thickness. 
Their  multiple  equilibria  included  three  different  anti-symmetric  equilibria  which  they 
called  type  Al,  A2,  and  A3,  as  well  as  two  pairs  of  non-symmetric  equilibria  which  they 
called  type  N1  and  N2. 

For  the  entire  region  of  parameter  space  explored  in  the  present  study,  only  one  anti¬ 
symmetric  equilibrium  was  found.  In  the  terminology  of  Cessi  and  lerley  this  equilibrium 
solution  is  said  to  be  of  type  A.  Other  regions  of  parameter  space  have  alternate  equilibria 
which  come  into  existence  via  bifurcations  of  this  anti-symmetric  equilibrium.  These 
equilibria  are  non-symmetric  and  come  in  pairs.  Three  types  of  non-symmetric  equilibria 
bifurcate  from  the  anti-symmetric  equilibrium,  those  of  type  N1  and  N2  as  well  as  a  third 
not  found  by  Cessi  and  lerley,  which  we  define  to  be  of  type  N3.  Finally,  the  solution 
branch  of  type  N1  undergoes  a  fold  catastrophe,  whereby  two  additional  non-symmetric 
equilibria  come  into  existence.  We  call  these  type  Nli  and  NI2.  Each  of  these  solution 
types  are  discussed  in  the  following  subsections. 

2.3.2  Anti-Symmetric  Equilibria  (type  A) 

In  Figure  2.1,  a  parameter  chart  indicates  the  bifurcations  of  the  anti-symmetric  equi¬ 
libria,  discussed  below.  Typical  anti-symmetric  stream-function  and  potential  vorticity 
fields  are  shown  in  Figure  2.2.  This  sequence  of  anti-symmetric  equilibria  is  taken  along 
the  left  most  side  of  the  parameter  chart  shown  in  Figure  2.1.  Apart  from  being  anti¬ 
symmetric,  these  equilibrium  solutions  are  characterized  by  the  formation  of  inertial 
recirculation  cells  flanking  the  southern  and  northern  edge  of  the  inter-gyre  boundary. 
It  is  important  to  mention  that  the  formation  of  closed  recirculation  cells  trapped  near 
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the  western  wall  does  not  occur  when  Ss  is  bigger  than  6h.  This  point  will  be  further 
discussed  in  Section  2.3.7.  The  remainder  of  this  section  is  restricted  to  the  case  Sh  >  Ss. 

For  moderate  values  of  5/,  the  recirculation  cells  are  small  and  trapped  near  the 
western  wall  (Figure  2.2,  Sj  =  0.03).  In  the  region  of  parameter  space  near  the  first 
pitchfork  bifurcation  (lower  part  of  the  lower  dashed  lobe  in  Figure  2.1),  both  viscous 
effects  and  advection  of  planetary  and  relative  vorticity  are  important  near  the  western 
boundary-layer  and  in  the  region  near  the  recirculation  gyres.  Note  how  the  q  contours 
run  almost  parallel  to  the  western  wall,  and  form  a  sharp  front  where  the  latitude  of 
zero  wind-stress  curl  intersects  the  western  boundary.  Away  from  the  western  boundary 
and  recirculation  cells,  the  flow  is  essentially  in  Sverdrup  balance.  In  this  region,  the  q 
contours  are  now  parallel  to  the  latitude  lines.  For  increasing  values  of  S[,  the  recircu¬ 
lations  cells  expand  in  size,  in  both  the  zonal  and  meridional  directions.  For  parameter 
values  close  to  the  second  pitchfork  bifurcation  (upper  part  of  the  lower  dashed  lobe  in 
Figure  2.1),  the  recirculation  gyres  have  expanded  across  the  basin  to  the  eastern  wall 
(Figure  2.2,  S[  =  0.06).  Further  increasing  <5/  leads  to  basin-filling  inertial  gyres  (Fig¬ 
ure  2.2,  Si  =  0.12).  Note  how  the  q  or  contours  are  parallel  to  the  stream  lines.  The 
Sverdrup  balance  is  then  no  longer  valid  anywhere  in  the  basin.  The  dominant  balance 
is  quasi-inertial,  i.e.  V^V')  ~  0-  Bottom  friction  which  is  not  important  in  any  local 
balance  is  nevertheless  important  in  maintaining  the  gyre  integrated  vorticity  balance. 

The  lower  dashed  lobe  in  Figure  2.1  is  the  location  of  bifurcation  points  at  which  non- 
symmetric  equilibria  of  type  N1  and  N2  bifurcate  from  the  anti-symmetric  equilibrium. 
The  lower  part  of  the  dashed  curve  gives  rise  to  steady  states,  (also  known  as  fixed 
points),  of  type  N1  and  the  upper  part  gives  rise  to  fixed  points  of  type  N2.  The  non- 
symmetric  equilibrium  of  type  N1  and  N2  will  be  described  in  Sections  2.3.3  and  2.3.4 
respectively.  The  upper  dashed  lobe  inside  the  dotted  lobe  in  Figure  2.1  is  the  location  of 
bifurcation  points  leading  to  non-symmetric  equilibria  of  type  N3.  Equilibria  of  type  N3 
are  described  in  Section  2.3.5.  Figure  2.5  which  is  fully  discussed  in  Section  2.3.3,  shows 
a  bifurcation  plot  of  ma.x{ip)  4-  min(^)  (a  measure  of  the  asymmetry  of  the  solution) 
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Figure  2.1:  Parameter  chart  for  5h  =  0.04,  showing  the  bifurcations  of  the  anti-symmetric 
equilibrium  (A).  Dashed  lines  indicate  pitchfork  bifurcations  and  dotted  lines  indicate 
Hopf  bifurcations.  The  anti-symmetric  equilibrium  exist  for  the  full  range  of  parameters 
shown  in  the  chart.  It  is  unstable  in  the  shaded  region,  and  stable  in  the  unshaded  region. 
The  lower  dashed  lobe,  indicates  where  non-symmetric  equilibria  of  type  N1  (lower  part 
of  lower  dashed  lobe),  and  of  type  N2  (upper  part  of  lower  dashed  lobe)  bifurcate  from  the 
anti-symmetric  equilibrium  via  a  symmetry  breaking  pitchfork  bifurcation.  The  upper 
dashed  lobe  indicates  where  non-symmetric  equilibria  of  type  N3  bifurcate  from  the  anti¬ 
symmetric  branch  via  a  pitchfork  bifurcation. 
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Figure  2.2:  Contour  plots  of  ip  (top  row)  and  q  (bottom  row)  for  the  anti-symmetric 
equilibrium  (type  A)  with  Ss  =  0.01,  Sh  =  0.04,  and  5j  as  indicated.  The  dashed  lines 
indicate  the  negative  contours  and  the  solid  lines  indicate  the  positive  contours.  The 
thick  solid  line  indicates  the  zero  contour.  The  contour  interval  is  0.03  for  ip  and  0.02  for 


versus  5/  for  a  value  of  Ss  =  0.01,  which  is  along  the  left  most  side  of  the  parameter  chart 
shown  in  Figure  2.1.  It  shows  the  four  symmetry  breaking  picthfork  bifurcations  that 
occur  as  one  moves  upwards  and  across  the  dashed  lobes  in  Figure  2.1. 

The  anti-symmetric  equilibrium  is  unstable  for  parameter  values  in  the  shaded  regions 
of  Figure  2.1,  and  stable  in  the  unshaded  region.  The  marginal  stability  curve  separating 
stable  from  unstable  regions  of  parameter  space  is  composed  of  a  curve  along  which  Hopf 
bifurcations  occur  (dotted  line)  and  a  curve  along  which  pitchfork  bifurcations  occur 
(lower  dashed  lobe).  Hopf  bifurcations  lead  to  self  sustained  oscillations  of  the  flow  field, 
and  pitchfork  bifurcations  are  at  the  origin  of  symmetry  breaking  multiple  equilibria. 

2.3.3  Non-symmetric  Equilibrium  (type  Nl) 

In  Figure  2.3,  a  parameter  chart  indicates  the  bifurcations  associated  with  the  non- 
symmetric  equilibrium  state  of  type  Nl.  Typical  non-symmetric  stream-function  and 
potential  vorticity  fields  of  type  Nl  are  shown  in  Figure  2.4,  along  with  equilibria  of  type 
A  and  N2  for  comparison.  As  mentioned  above,  the  non-symmetric  states  of  type  Nl 
come  into  existence  via  a  symmetry  breaking  pitchfork  bifurcation  of  the  anti-symmetric 
state.  The  pitchfork  bifurcation  points  that  mark  the  emergence  of  the  type  Nl  equilibria 
are  located  on  the  lower  dashed  curve  in  Figure  2.3.  The  non-symmetric  equilibria  exist 
only  in  the  region  bounded  by  the  lower  dashed  curve,  and  the  upper  solid  curve.  Within 
this  region,  it  is  stable  in  the  unshaded  area,  and  unstable  in  the  shaded  area.  Within 
the  wedge  shaped  region  bounded  by  the  saddle  node  bifurcations  emanating  from  the 
cusp  point  labeled  CPI  there  are  two  additional  equilibria  defined  to  be  of  type  Nli 
and  NI2.  The  distinction  between  equilibria  of  types  Nl,  Nli,  and  NI2  is  essentially 
in  the  geometrical  arrangement  of  the  multiple  circulation  cells  within  the  basin.  A  full 
discussion  of  equilibria  of  type  Nli,  and  NI2  will  be  given  in  Section  2.3.6.  Along  the 
solid  curve  labeled  SNl,  the  Nl  branch  experiences  a  saddle-node  bifurcation  and  non- 
symmetric  equilibria  of  type  Nli  come  into  existence.  The  equilibria  of  type  Nli  also 
experience  a  saddle-node  bifurcation  along  the  lower  solid  curve  labeled  SN2,  and  non- 
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symmetric  equilibria  of  type  NI2  come  into  existence.  Refering  again  to  Figure  2.5,  one 
can  see  the  saddle-node  bifurcations  as  one  moves  across  the  solid  curves  in  Figure  2.3. 

In  a  neighborhood  of  the  pitchfork  bifurcation  point  the  equilibria  are  characterized 
by  moderate  recirculation  gyres  with  essentially  the  same  vorticity  balance  as  the  anti¬ 
symmetric  state,  provided  that  6n  >  5s,  see  discussion  in  Section  2.3.7.  As  5[  is  increased, 
and  one  moves  away  from  the  bifurcation  point,  one  of  the  recirculation  cells  becomes 
stronger.  The  weaker  of  the  two  recirculation  cells  crosses  the  line  of  zero  wind-stress 
curl  as  it  is  pulled  by  its  more  intense  counterpart.  As  5s  is  decreased,  and  5/  is  kept 
near  the  center  of  the  shaded  region  in  Figure  2.3,  the  flow  field  becomes  progressively 
more  asymmetric,  with  large  meanders  crossing  the  latitude  of  zero  wind-stress  curl. 
These  meanders  permit  vorticity  to  be  fluxed  across  the  latitude  of  zero  wind-stress 
curl,  so  that  the  integrated  vorticity  balance  need  not  be  achieved  entirely  by  bottom 
friction.  This  allows  equilibria  of  type  N1  to  remain  of  much  weaker  intensity  than  its 
anti-symmetric  counterpart  of  type  A,  as  5/  is  increased.  Equilibria  of  type  A,  Nl,  and 
N2  can  be  compared  in  Figure  2.4.  One  should  note  in  particular  how  the  non-symmetric 
equilibrium  of  type  Nl  are  somewhat  intensified  in  the  western  part  of  the  basin,  and 
much  weaker  than  those  of  type  A  or  N2.  A  further,  and  more  dramatic  example  of 
the  different  vorticity  balance  achieved  by  non-symmetric  equilibria  of  type  Nl,  can  be 
observed  in  Figure  2.6,  by  comparing  equilibria  of  type  A,  Nl,  and  N2.  The  multiple 
equilibria  in  this  figure  were  computed  for  the  same  parameter  values  as  those  used 
in  the  time-dependent  simulation  of  Marshall  (1984),  (5s  —  0.001,  5h  =  0.0313,  and 
5[  =  0.0316).  The  stream-function  field  labeled  PBAR  in  Figure  2.6  is  the  time-average 
of  a  simulation  of  the  time-dependent  flow  field.  The  equilibria  of  type  Nl  are  unstable  in 
this  region  of  parameter  space.  The  time-dependent  flow  field  evolved  in  a  complicated 
way  with  many  strong  eddies  forming.  Nevertheless,  the  integrated  vorticity  balance  over 
the  region  of  negative  wind-stress  curl  is  similar  to  that  of  time-averaged  flow,  provided 
the  role  of  the  eddy  induced  vorticity  flux  in  the  time-dependent  case  is  replaced  by  the 
flux  of  vorticity  across  the  line  of  zero  wind-stress  curl  by  stationary  meanders  in  the 
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Figure  2.3:  Parameter  chart  for  Sjj  =  0.04,  showing  the  bifurcations  of  the  non-symmetric 
equilibria  of  type  Nl.  Solid  lines  denote  saddle-node  bifurcations,  dashed  lines  indicate 
pitchfork  bifurcations,  and  dotted  lines  indicate  Hopf  bifurcations.  Equilibria  of  type  Nl 
can  be  traced  continuously  for  increasing  values  of  5/  from  the  dashed  curve  (where  they 
bifurcate  from  the  anti-symmetric  equilibrium),  up  to  the  saddle-node  curve  labeled  SN 
or  SNl  depending  on  whether  6s  is  to  the  right  or  left  of  the  cusp  point  labeled  CPI. 
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Figure  2.5:  Bifurcation  plot  of  V’max  +  '^min  -vs-  <5/,  indicating  the  emergence  of  non- 
symmetric  equilibria  N1  and  N2  and  N3  via  pitchfork  bifurcations  at  Sj  =  5pi,  5p2,  and 
5p3.  A  fourth  pitchfork  bifurcation  point  at  6i  =  5p^  marks  the  disappearance  of  the 
non-symmetric  equilibria  of  type  N3.  Anti-symmetric  equilibria  lie  on  the  horizontal  line 
in  the  center  of  the  figure.  Saddle-node  bifurcation  points  at  8[  =  ^51,  ^52,  and  ^53, 
mark  the  merging  of  equilibria  of  type  N1  with  Nli,  Nli  with  NI2,  and  of  type  NI2  with 
N2.  The  solid  curves  indicate  equilibria  which  are  stable  and  the  dashed  lines  indicate 
equilibria  which  are  unstable.  {83  =  0.01,  8h  =  0.04). 
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Figure  2.6:  Contour  plots  of  V’  for  5i  =  0.0316,  5s  =  0.001  and  5h  =  0.0313.  Steady  states 
of  type  A,  Nl,  and  N2  are  presented  along  with  the  time-averages  stream-function  field 
(PBAR).  The  solid  lines  denote  positive  contours  and  the  dashed  lines  denote  negative 
contours.  Contour  intervals  are  as  indicated. 

steady  case.  In  the  steady  state  solution,  bottom  friction  removed  only  4  percent  of  the 
vorticity  input  by  the  wind  and  in  the  time-mean  state,  bottom  friction  removed  about  5 
percent  of  the  vorticity  input  by  the  wind.  Furthermore,  the  pair  of  fixed  points  of  type 
Nl  are  qualitatively  similar  to  a  typical  instantaneous  flow  field  of  the  time-dependent 
simulation. 

•  For  5s  >  0.0161  the  non-symmetric  equilibrium  of  type  Nl  can  be  traced  continu¬ 
ously  up  to  a  saddle-node  bifurcation  curve  (labeled  SN  in  Figure  2.3),  where  it  merges 
with  a  fixed  point  of  type  N2  and  ceases  to  exist.  At  5s  =  0.0161,  and  5i  =  0.090,  a 
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fold  catastrophe  occurs  (labeled  CPI  in  Figure  2.3),  which  leads  to  a  wedge  with  two 
additional  equilibrium  states  called  type  Nli  and  NI2.  These  additional  equilibria  will 
be  described  in  Section  2.3.6.  Furthermore,  for  Ss  <  0.0161,  i.e.  for  values  of  63  to  the 
left  CPI,  it  is  NI2  that  merges  with  N2  at  SN. 

2.3.4  Non-Symmetric  Equilibrium  (type  N2) 

In  Figure  2.7,  a  parameter  chart  indicates  the  bifurcations  associated  with  the  non- 
symmetric  equilibrium  state  of  type  N2.  Typical  non-symmetric  stream-function  and 
potential  vorticity  fields  of  type  N2  were  shown  in  Figure  2.4.  As  mentioned  in  Sec¬ 
tion  2.3.1,  non-symmetric  equilibria  of  type  N2  come  into  existence  via  a  pitchfork  bi¬ 
furcation  of  the  anti-symmetric  equilibrium.  This  pitchfork  bifurcation  has  its  origin  at 
the  cusp  point  labeled  CP  in  Figure  2.7.  Equilibria  of  type  N2  are  always  unstable.  This 
branch  can  be  traced  continuously  up  to  a  saddle-node  bifurcation  point  lying  on  the  solid 
curve,  where  it  merges  with  the  equilibrium  of  type  N1  and  disappears.  The  pitchfork 
bifurcation  and  the  saddle-node  bifurcation  can  be  seen  in  Figure  2.5  for  Ss  =  0.01. 

Near  the  cusp  point  in  Figure  2.7,  there  is  little  difference  between  equilibria  of  type 
N1  and  N2.  In  this  region  of  parameter  space  a  three  term  balance  is  achieved  between 
the  advection  of  planetary  vorticity,  advection  of  relative  vorticity,  and  lateral  diffusion. 
However,  the  vorticity  balance  integrated  over  the  region  of  either  positive  or  negative 
wind-stress  curl  is  achieved  primarily  by  bottom  friction.  As  Ss  is  decreased,  and  Sj 
is  kept  in  the  center  of  the  shaded  region  in  Figure  2.7,  the  equilibria  of  type  N1  and 
N2  separate  themselves  through  different  balances.  The  non-symmetric  equilibrium  of 
type  N2  is  characterized  by  inertial  gyres  V^V’)  ~  0),  that  extend  across  the  basin. 
Again,  the  bottom  friction  term  is  unimportant  in  any  local  balance,  but  is  nonetheless 
crucial  in  maintaining  the  global  vorticity  balance.  This  should  be  contrasted  with  non- 
symmetric  equilibria  of  type  Nl,  where  viscous  effects  are  important  locally  in  the  region 
near  the  western  wall,  and  where  the  inter-gyre  boundary  meets  the  western  wall. 
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Figure  2.7:  Parameter  chart  for  =  0.04,  showing  the  bifurcations  of  the  non-symmetric 
equilibrium  (N2).  Solid  lines  denote  saddle-node  bifurcations  and  dashed  lines  indicate 
pitchfork  bifurcations.  Equilibria  of  type  N2  exist  only  in  the  shaded  region  of  param¬ 
eter  space  where  it  is  also  always  unstable.  Non-symmetric  equilibria  of  type  N2  come 
into  existence  via  a  pitchfork  bifurcation  of  the  anti-symmetric  flow,  and  can  be  traced 
continuously  as  Si  is  increased  up  to  the  solid  curve  labeled  SN  where  it  experiences  a 
saddle-node  bifurcation,  merges  with  equilibria  of  type  Nl,  and  ceases  to  exist. 
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2.3.5  Non-symmetric  Equilibrium  (type  N3) 

Non-symmetric  equilibria  of  type  N3  exist  inside  the  region  bounded  by  the  upper  dashed 
lobe  of  Figure  2.1.  Typical  xp  and  q  fields  for  the  non-symmetric  equilibria  of  type  N3  and, 
for  comparison,  anti-symmetric  equilibria  are  presented  in  Figure  2.8.  Non-symmetric 
equilibria  of  type  N3  come  into  existence  via  a  pitchfork  bifurcation  of  the  anti-symmetric 
equilibrium.  The  equilibria  is  characterized  by  basin-filling  inertial  gyres  with  the  inter¬ 
gyre  boundary  shifted  either  north  or  south,  and  tilted  in  the  east-west  direction.  Both 
equilibria  of  type  A  and  N3  are  characterized  by  a  quasi-inertial  balance,  the  difference 
in  the  solution  is  restricted  to  their  different  symmetry  properties. 

2.3.6  Non-Symmetric  Equilibria  (type  Nli  and  NI2) 

Non-symmetric  equilibria  of  type  Nli  and  NI2  come  into  existence  via  a  fold  catastrophe 
of  the  solution  branch  of  type  Nl.  In  Figure  2.3,  the  region  of  parameter  space  where  this 
fold  takes  place  is  denoted  by  the  wedged  shaped  region  emanating  from  the  cusp  point 
labeled  CPI.  Figure  2.9  shows  typical  contour  plots  of  ijj  and  q.  These  equilibria  are 
characterized  by  strongly  non-symmetric  flow  fields,  with  multiple  circulation  cells.  The 
essential  difference  between  equilibria  of  type  Nl,  Nli,  and  NI2,  is  characterized  by  the 
geometrical  arrangement  of  circulation  cells  within  the  basin.  The  integrated  vorticity 
balance  for  the  region  of  either  positive  or  negative  wind-stress  curl  is  achieved  primarily 
by  stationary  meanders  transporting  vorticity  across  the  latitude  of  zero  wind-stress  curl. 

2.3.7  Dependence  on  Lateral  DifFusivity  Parameter,  Sh 

In  this  section,  the  role  of  the  lateral  diffusion  parameter,  in  modifying  the  model’s 
fixed  points  is  explored.  It  is  interesting  to  compare  two  sequences  of  anti-symmetric 
equilibrium  solutions  with  alternate  orderings  of  lateral  friction  and  bottom  friction  layer 
thicknesses,  i.e.  Sg  >  Ss  and  63  >  Sh-  In  this  section,  computations  with  6h  =  0.01 
were  performed  with  TVj.  =  33,  and  Ny  =  129,  grid  points  on  a  grid  stretched  in  the  y 
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Figure  2.8;  Typical  and  q  fields  for  the  branches  Al  and  N3.  6j  =  0.1127,  Ss  =  0.01 
and  Sff  =  0.04.  The  contour  interval  is  0.4  for  ip  and  0.3  for  q.  The  negative  contours 
are  dashed,  and  the  zero  contour  is  the  thick  one.  The  non-symmetric  branch  N3  has  its 
inter-gyre  boundary  shifted  northward  and  tilted  from  west  to  east. 


Figure  2.9:  Typical  xp  and  q  fields  for  the  branches  Nl,  Nli,  and  NI2,  for  61  =  0.1, 
Ss  =  0.02  and  ~  0.04.  The  contour  interval  is  0.2  for  xp  and  0.2  for  q.  The  negative 
contours  are  dashed,  and  the  zero  contour  is  the  thick  one.  The  difference  between 
the  three  equilibrium  branches  is  essentially  restricted  to  the  spatial  arrangement  of  the 
multiple  closed  circulation  cells. 
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direction  with  with  rriy  =  0.25.  See  Appendix  B  for  a  definition  of  nix  and  niy.  The  grid 
spacing  in  the  x  direction  was  uniform  dx  =  0.03,  and  in  the  y  direction  the  grid  spacing 
varied  from  a  maximum  of  dy  =  0.03  near  the  northern  and  southern  basin  walls,  and 
a  minimum  of  dy  =  0.0031,  near  the  line  of  zero  wind-stress  curl.  The  higher  resolution 
and  stretched  grid  are  necessary  to  resolve  the  sharp  gradient  in  the  relative  vorticity 
field  that  forms  where  the  inter-gyre  boundary  meets  the  western  wall. 

Figure  2.2  shows  a  sequence  of  anti-symmetric  equilibria  for  increasing  values  of  5/ 
with  5h  =  0.04  >  Ss  =  0.01  fixed.  For  this  ordering  of  the  friction  parameters,  the 
sequence  of  stream-function  fields  is  similar  to  that  computed  by  lerley  and  Sheremet 
(1995)  for  the  single-gyre  case  with  Munk  type  lateral  diffusion.  Closed  recirculation  cells 
form  near  the  western  wall  where  the  counter  rotating  gyres  meet.  The  cells  expand  in 
size  and  strength  as  5i  is  increased,  eventually  filling  the  entire  basin. 

Contrast  this  with  Figure  2.10  which  shows  a  similar  sequence  but  with  the  alternate 
ordering  in  the  thickness  of  the  lateral  and  bottom  friction  layers,  Sh  =  0.01  <  Ss  =  0.04. 
In  this  case,  the  sequence  of  equilibria  is  similar  to  the  sequence  of  solutions  computed 
by  Veronis  (1966).  The  jet  separating  the  northern  gyre  from  the  southern  gyre  increases 
in  strength  and  penetrates  progressively  deeper  into  the  interior  as  Si  increases.  The 
circulation  pattern  does  not  develop  closed  recirculation  cells  trapped  near  the  western 
wall.  The  limit  <5/  — >  oo,  has  basin-filling  gyres  similar  to  those  for  the  case  with  Sh  >  Ss, 
in  the  sense  that  they  both  have  q  versus  ip  scatter  plots  with  negative  slope.  For  the 
case  of  Ss  >  Sh,  there  are  however  intermediate  values  of  Sj  for  which  the  scatter  plot  of 
q  vs.  'll)  has  positive  slope. 

Finally,  in  comparing  the  single-gyre  calculations  of  Veronis  (1966),  to  the  double¬ 
gyre  calculations  presented  here,  one  should  highlight  the  importance  of  the  boundary 
condition  experienced  by  the  northern  flank  of  the  sub-tropical  gyre.  For  a  single-gyre 
calculation,  there  can  be  no  flux  of  vorticity  across  the  northern  wall,  which  is  not  the 
case  for  the  double-gyre  calculations.  This  difference  in  boundary  condition  allows  the 
northern  jet  to  penetrate  clear  across  the  basin  for  either  ordering  of  5//  and  Ss-  Tight 
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Figure  2.10;  Contour  plots  of  ^  (top  row)  and  q  (bottom  row)  for  the  anti-symmetric 
equilibrium  (type  A)  with  8s  =  0.04  and  8h  =  0.01  and  5i  as  indicated.  The  dashed  lines 
indicate  the  negative  contours  and  the  solid  lines  indicate  the  positive  contours.  The 
thick  solid  line  indicates  the  zero  contour.  C.I.  =  0.02  for  ip  and  C.I.  =  0.02  for  q. 
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Figure  2.11:  Contour  plots  of  V'  and  q  for  the  single-gyre  case.  The  no-flux  of  vorticity 
across  the  northern  wall  allows  the  northern  jet  to  penetrate  straight  across  the  basin. 
This  should  be  compared  with  the  double-gyre  anti-symmetric  equilibria  computed  for 
the  same  parameters  (Figure  2.2,  5[  —  0.03).  5s  =  0.01  and  6h  =  0.04  and  5i  =  0.03. 
C.I.  =  0.02  for  V’  and  C.I.  =  0.02  for  q. 


recirculation  cells  trapped  in  the  north  western  corner  of  the  sub-tropical  gyre  can  form 
only  if  5h  >  5$  and  anomalously  high  potential  vorticity  is  allowed  to  diffuse  from  the 
sub-polar  gyre  into  the  sub-tropical  gyre.  To  demonstrate  this,  Figure  2.11  shows  contour 
plots  of  the  q  and  xjj  field  for  a  single-gyre  calculation  with  Si  =  0.03,  Sh  =  0.04  and 
Ss  ~  0.01.  This  single-gyre  calculation  should  be  compared  with  the  double-gyre  steady 
state  equilibrium  for  the  same  parameter  values  (Figure  2.2).  Incidentally,  the  single-gyre 
model  with  5h  =  0.04  and  Ss  =  0.01  does  not  exhibit  a  saddle  node  bifurcation  as  Si  is 
increased.  This  is  different  from  the  results  of  lerley  and  Sheremet  (1995).  It  is  not  clear 
whether  their  bifurcation  occurs  for  different  Ss  and  Sh,  or  if  the  diflferent  dissipation 
operator  and  boundary  condition  prevent  it  from  occurring  in  the  present  model. 

At  least  near  the  region  of  parameter  space  explored,  the  bifurcation  structure  leading 
to  non-symmetric  equilibria  of  type  N1  and  N2  does  not  depend  qualitatively  on  the 
value  of  Sh-  Figure  2.12  shows  a  plot  of  max{il))  4-  min{ip)  (a  measure  of  asymmetry) 
vs.  Si,  for  Ss  =  0.04,  and  Sh  set  at  either  =  0.04  (solid  line),  or  0.01  (dashed  line). 
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Figure  2.12:  Bifurcation  plot  of  maxV^  +  minV’  -vs-  indicating  the  emergence  of  non- 
symmetric  equilibria  of  type  N1  and  N2,  via  pitchfork  bifurcations  for  two  different  values 
of  J//,  {5h  =  0.04  solid  curve,  5n  =  0.01  dotted  curve),  5$  —  0.04,  and  their  disappearance 
at  saddle-node  bifurcations. 


The  plot  shows  the  emergence  of  non-symmetric  equilibria  of  type  Nl  and  N2  as  S[  is 
increased.  For  =  0.01  the  pitchfork  bifurcations  occur  for  slightly  smaller  values 
of  Sj  than  they  do  for  6h  =  0.04.  Also  for  the  smaller  value  of  6h,  the  saddle-node 
bifurcation,  where  the  equilibria  of  type  Nl  and  N2  merge,  happens  at  a  larger  value  of 
6[.  The  top  row  of  Figure  2.13  shows  contour  plots  of  the  ip  field,  and  the  bottom  row 
shows  contour  plots  of  the  q  field  as  one  moves  around  the  nose  from  point  marked  1 
through  5  in  Figure  2.12.  The  sequence  begins  with  the  anti-symmetric  equilibria  at  the 
location  where  Nl  bifurcates  from  the  anti-symmetric  solution  (location  1),  continues 
with  a  typical  non-symmetric  equilibrium  of  type  Nl  (location  2),  the  non-symmetric 
equilibrium  at  the  nose  point  NP  where  equilibria  of  type  Nl  and  N2  merge  (location  3), 
a  non-symmetric  equilibrium  of  type  N2  (location  4),  and  ends  with  the  anti-symmetric 
equilibrium  at  the  point  where  equilibria  of  type  N2  bifurcates  from  the  anti-symmetric 
branch  (location  5).  Note  that  the  tight  recirculation  cell  that  form  near  the  western 
wall  for  the  case  5h  >  Ss  are  not  present  when  Ss  >  Sh- 

2.4  Overview  of  Bifurcation  Structure 

It  is  useful  to  summarize  the  bifurcation  structures  described  in  the  previous  sections.  In 
Figure  2.5,  the  difference  in  the  extreme  values  of  the  magnitude  of  the  stream-function 
field  in  the  sub-polar  and  sub-tropical  gyres  {ipmax  +  i’min)  is  plotted  as  a  function  of 
Si-  This  bifurcation  plot  is  taken  along  the  left  most  edge  of  the  parameter  charts 
in  Figures  2.1,  2.3,  and  2.7,  and  thus  includes  all  the  bifurcations  described  in  the 
previous  sections.  It  illustrates  the  successive  symmetry  breaking  bifurcations  of  the  anti¬ 
symmetric  branch.  The  anti-symmetric  branch  is  marked  by  the  horizontal  line  in  the 
center  of  the  plot.  A  pair  of  non-symmetric  equilibria  of  type  Nl  emerges  as  5/  is  increased 
past  the  first  pitchfork  bifurcation  point  at  5pi.  The  pair  of  non-symmetric  equilibria 
of  type  N2  emerges  as  is  increased  past  the  second  pitchfork  bifurcation  at  Sp2.  The 
symmetry  breaking  equilibria  of  type  N3  emerges  as  S[  is  increased  past  the  pitchfork 
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Figure  2.13:  Contour  plots  of  'ip  (top  row)  and  q  (bottom  row)  around  the  nose  point 
where  the  branches  N1  and  N2  merge.  The  solutions  correspond  to  values  of  =0.07017, 
0.08041,  0.1429,  0.1245,  0.1087  and  Ss  =0.04.  The  first  solution  labeled  Pf  is  close  to  the 
first  pitchfork  bifurcation.  The  second  solution  labeled  N1  is  a  typical  solution  of  type 
Nl.  The  third  solution  is  close  to  the  saddle- node  bifurcation  where  N1  and  N2  merge. 
The  fourth  solution  is  a  typical  solution  of  type  N2,  and  the  last  solution  is  close  to  the 
second  pitchfork  bifurcation.  The  dashed  lines  indicate  the  negative  contours  and  the 
solid  lines  indicate  the  positive  contours.  The  thick  solid  line  indicates  the  zero  contour. 
C.I.  =  0.02  for  ip  and  C.I.  =  0.02. 
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bifurcation  point  at  ^P3.  As  5/  is  increased  past  the  fourth  pitchfork  bifurcation  point 
at  <5p4,  the  pair  of  non-symmetric  equilibria  of  type  N3  merges  with  the  anti-symmetric 
branch  and  disappear.  The  plot  also  shows  the  saddle-node  bifurcation  point  where  the 
equilibria  of  type  Nli  and  NI2  are  created,  {61  =  5si  and  832),  as  well  as  the  saddle-node 
bifurcation  point  (5i  =  ^53)  where  the  equilibria  of  type  NI2  and  N2  merge. 

Figure  2.14  is  a  composite  figure  displaying  an  overlay  of  all  the  bifurcation  curves. 
In  each  region  of  parameter  space  a  pair  of  numerals  indicates  the  total  number  of  steady 
equilibria  coexisting  for  the  same  parameter  values  and  the  number  of  those  which  are 
stable.  For  example  (5,1)  would  imply  that  there  are  5  equilibrium  states  one  of  which 
is  stable.  Regions  with  multiple  stable  equilibria  are  limited  to  those  marked  by  (3,2) 
where  the  pair  of  equilibria  of  type  N1  are  stable,  and  the  anti-symmetric  equilibrium  is 
unstable,  as  well  as  to  the  region  denoted  by  (5,3)  where  the  stable  equilibria  are  of  type 
A  and  Nl,  and  the  unstable  equilibria  are  of  type  N2.  Regions  with  up  to  11  unstable 
equilibria  are  identified  in  the  figure. 

2.5  Discussion 

Using  a  continuation  algorithm  for  finding  both  steady  state  solutions,  and  their  cor¬ 
responding  least  stable  eigenmode  (or  one  unstable  eigenmode  if  it  exists),  we  mapped 
the  bifurcation  structure  of  the  steady  state  solutions  of  a  barotropic  wind-driven  ocean 
model  as  a  function  of  the  two  nondimensional  parameters  —  the  inertial  layer  thickness, 
S[,  and  the  Stommel  layer  thickness,  63-  One  of  the  goals  in  carrying  out  these  calcula¬ 
tions  was  to  contribute  to  the  broader  objective  of  mapping  the  states  of  ocean  models 
with  a  varying  complexity  of  sub-grid-scale  parameterizations,  and  boundary  conditions. 
The  model  we  used  had  bottom  friction  and  lateral  diffusion  with  super-slip  boundary 
conditions.  In  this  sense  this  study  is  a  continuation  of  the  work  of  Cessi  and  lerley 
(1995). 
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Figure  2.14:  Parameter  chart  for  5h  =  0.04  showing  an  overlay  of  the  bifurcations  for 
all  the  branches  found.  Solid  lines  denote  saddle-node  bifurcations,  dashed  lines  indicate 
pitchfork  bifurcations,  and  dotted  lines  indicate  Hopf  bifurcations. 


2.5.1  Inertial  runaway 

One  solution  branch  (type  A)  can  be  traced  continuously  from  the  linear  regime  to  the 
highly  inertial  regime.  Like  the  forcing  function,  this  branch  is  anti-symmetric  about  the 
mid-basin  latitude.  For  sufficiently  strong  forcing  or  sufficiently  weak  bottom  friction,  the 
anti-symmetric  solution  tended  towards  a  highly  inertial  circulation  with  transports  far 
in  excess  of  those  predicted  by  Sverdrup  balance.  Apart  from  having  unrealistically  large 
transports,  the  inertial  runaway  solution  shows  no  westward  intensification.  Furthermore, 
for  all  values  of  5s  explored,  increasing  Sj  eventually  leads  to  a  region  where  the  anti¬ 
symmetric  steady  state  is  stable,  and  apparently  unique.  A  limited  number  of  time- 
dependent  calculations  with  parameters  in  the  region  labeled  (1,1)  in  Figure  2.14  all 
converged  to  the  anti-symmetric  fixed  point,  regardless  of  the  initial  conditions  used, 
suggesting  that  this  fixed  point  is  a  global  attractor  at  sufficiently  large  values  of  (5/. 

The  stability  of  the  anti-symmetric  solution  as  5i  oo,  should  be  contrasted  with  the 
stability  results  of  Cessi  and  lerley  (1995),  for  the  nonlinear  Munk  model  with  free-slip 
boundary  conditions  at  the  eastern  and  western  walls.  They  found  that  for  <5/  oo  and 
5m  ‘C  1,  where  5m  is  the  nondimensional  Munk  layer  thickness,  the  only  equilibrium 
is  unstable,  and  has  a  single  unstable  eigenmode.  The  difference  in  stability  between 
Cessi  and  lerley’s  model  and  the  one  considered  in  this  study  can  be  attributed  to  the 
bottom  friction  term  as  opposed  to  the  choice  of  lateral  diffusion  operator  and  boundary 
conditions.  To  verify  this,  some  of  the  calculations  of  Cessi  and  lerley  were  repeated  with 
and  without  bottom  friction.  The  calculations  were  carried  out  to  values  of  5[  as  high 
as  300,  with  5s  =  0.  The  growth-rate  of  the  unstable  eigenmode  remained  positive  but 
decreased  monotonically.  The  computations  were  then  repeated  with  the  addition  of  a 
bottom  friction  term,  with  a  finite  value  of  5s,  and  the  unstable  mode  became  stable  for 
a  sufficiently  large  value  of  5/. 

Another  difference  between  our  results  and  those  of  Cessi  and  lerley  (1995),  is  the 
non-existence  of  a  cusp  catastrophe  leading  to  multiple  anti-symmetric  equilibria.  Recall 
that  the  model  used  here  has  bottom  friction,  lateral  diffusion  in  the  form  of  a  biharmonic 
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operator  acting  on  the  vorticity,  and  super-slip  boundary  conditions  which  do  not  allow  a 
flux  of  vorticity  through  the  basin  walls.  We  have  investigated  whether  the  choice  between 
free-slip  and  super-slip  boundary  conditions  as  opposed  to  the  form  of  the  dissipation 
operator  could  be  responsible  for  the  existence  of  this  cusp  catastrophe.  Using  free-slip 
boundary  conditions  (C  =  V^C  =  0  on  dV)  we  deduced  the  existence  of  a  cusp  for  6s 

between  0.01  and  0.001.  On  the  other  hand,  for  the  super-slip  case,  no  cusp  leading 

\ 

to  multiple  anti-symmetric  equilibria  exists  for  values  of  bottom  friction  greater  than 
6s  >  0.001.  It  is  not  clear  if  the  cusp  exists  for  smaller  values  of  6s  and  6h,  or  if 
the  super-slip  boundary  conditions  truly  prevents  it  from  forming  for  any  value  of  the 
dissipation  parameters.  Using  a  single-gyre  model,  lerley  and  Sheremet  (1995),  found 
the  cusp  catastrophe  when  lateral  diffusion  and  free-slip  boundary  conditions  were  used. 
However,  they  did  not  find  the  cusp  catastrophe  when  bottom  friction  was  used  instead 
of  lateral  diffusion.  A  model  with  bottom  friction  alone  does  not  allow  vorticity  to  diflFuse 
through  the  basin  wall  and  is  similar,  in  this  sense,  to  a  model  with  lateral  diffusion  and 
a  super-slip  boundary  condition.  Perhaps  lateral  diffusion  of  vorticity  through  the  basin 
wall  is  an  essential  element  to  the  dynamical  balance  that  allows  for  the  existence  multiple 
anti-symmetric  equilibria.  Reducing  bottom  friction  and  holding  lateral  diffusion  fixed 
at  values  of  6ii  >  0.01  does  not  lead  to  a  fold  of  the  anti-symmetric  branch.  We  have 
not,  however,  eliminated  the  possibility  that  a  reduction  of  lateral  diffusion  would  give 
rise  to  a  fold  in  the  anti-symmetric  branch. 

2.5.2  Internal  Compensation 

As  already  pointed  out,  for  high  boundary-layer  Reynolds  number,  the  anti-symmetric 
solution  has  a  basin-filling  inertial  gyre  with  no  western  intensification  and  transports 
far  in  excess  of  those  predicted  from  Sverdrup  balance.  As  in  the  case  of  the  single-gyre 
solutions  found  by  lerley  and  Sheremet  (1995),  a  large  eddy  viscosity  is  needed  for  the 
model’s  anti-symmetric  solution  to  have  western  intensification  and  a  mass  transport 
comparable  to  that  observed  in  the  real  ocean.  lerley  and  Sheremet  (1995),  point  out 
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that  a  transfer  of  vorticity  between  adjacent  gyres  —  a  mechanism  which  is  precluded 
in  a  single-gyre  model  —  should  reduce  the  need  for  a  large  eddy  viscosity.  Indeed, 
in  idealized  double-gyre  models,  the  inter-gyre  transfer  of  vorticity  can  be  the  primary 
mechanism  by  which  the  wind-stress  curl  vorticity  input  is  balanced  within  each  half  of 
the  basin  as  shown  in  Harrison  and  Holland  (1981)  and  Marshall  (1984). 

For  double-gyre  models  that  have  no  net  input  of  vorticity  over  the  entire  basin,  there 
is  the  possibility  that  the  vorticity  balance  can  be  achieved  internally  without  the  need 
for  vorticity  to  be  fluxed  through  the  basin  walls  or  bottom.  Cessi  and  lerley  (1995), 
have  demonstrated  that  double-gyre  models  with  anti-symmetric  wind-stress  curl  profile 
admit  non-symmetric  solutions.  We  also  found  symmetry  breaking  pitchfork  bifurcations 
leading  to  non-symmetric  equilibria  (type  N1,N2,N3).  Because  of  their  non-symmetry 
with  respect  to  the  wind  forcing,  these  solutions  have  currents  which  transport  vortic¬ 
ity  across  the  latitude  of  zero  wind-stress  curl.  For  these  non-symmetric  solutions,  the 
input  of  vorticity  by  the  wind-stress  curl  in  each  half  of  the  basin  is  balanced  in  part 
by  the  export  of  vorticity  by  the  current.  This  additional  mechanism  for  removing  vor¬ 
ticity  lessens  the  burden  on  the  explicit  dissipation  for  balancing  the  vorticity  budget. 
Consequently  these  solutions  have  weaker  currents  than  the  anti-symmetric  solution. 
Figure  2.15  shows  a  plot  of  the  difference  between  the  maximum  and  minimum  of  the 
stream-function,  a  quantity  proportional  to  the  maximum  transport  in  the  basin,  as  a 
function  of  5/.  The  solid  curves  correspond  to  equilibria  that  are  stable  and  the  dashed 
lines  correspond  to  equilibria  that  are  unstable.  The  crosses  connected  by  the  dotted 
line  correspond  to  a  series  of  time-averaged  solutions.  The  initial  condition  for  these 
simulations  were  the  steady  state  of  type  N1  plus  some  random  noise  perturbation.  One 
can  see  that  the  transport  for  the  anti-symmetric  solution  (labeled  A),  is  always  higher 
than  that  of  the  non-symmetric  solutions  (labeled  Nl,  N2,  and  N3).  The  time-averaged 
solutions  also  show  reduced  transports  compared  to  the  anti-symmetric  steady  state  so¬ 
lutions.  In  fact,  the  time-averaged  transports  are  of  a  magnitude  comparable  to  those 
of  the  non-symmetric  branch  of  type  Nl.  It  appears  that  unstable  fixed  point  of  type 
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Figure  2.15;  Maximum  transport  for  steady  state  solution  and  time-averaged  solutions 
as  a  function  of  Sj.  Dotted  lines  with  indicate  time-averaged  solutions.  Solid  lines 
indicate  stable  steady  state  solutions  and  dashed  lines  indicate  unstable  steady  state 
solutions.  8s  =  0.01  and  8h  =  0.04. 

N1  are  close  to  a  more  complicated  attractor  which  includes  the  model  trajectories  for 
which  the  time-averaged  solutions  were  computed.  Whether  or  not  these  fixed  points  are 
part  of  this  attractor  remains  a  conjecture  since  we  have  not  attempted  to  compute  more 
complicated  orbits  like  limit  cylces  or  homoclinic  and  heteroclinic  orbits. 

The  non-symmetric  equilibria  only  exist  for  a  finite  range  of  parameters,  between 
the  first  pitchfork  bifurcation  point  labeled  PFl  and  the  last  turning  point  labeled  SN. 
Beyond  the  turning  point,  all  the  computed  model  trajectories  converged  to  the  anti¬ 
symmetric  inertial  runaway  solution.  Even  before  the  nose  point  is  reached,  some  model 
trajectories  asymptote  to  the  stable  anti-symmetric  branch.  In  the  regions  labeled  (5,1) 
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in  Figure  2.14,  it  is  the  initial  condition  that  determines  whether  the  model  trajectory 
fluctuates  in  some  complicated  manner  in  the  neighborhood  of  the  N1  fixed  points  or 
asymptotes  to  the  stable  equilibrium  of  type  A,  which  is  where,  as  already  mentioned, 
all  model  trajectories  converged  for  parameter  values  beyond  SN.  We  can  conclude  from 
these  results  that  inertial  runaway  is  unavoidable  as  the  boundary-layer  Reynolds  number, 
becomes  sufficiently  large. 

Inter-gyre  fluxes  of  vorticity  prevent  the  time-averaged  solution  from  developing  un¬ 
realistically  large  transports,  but  they  do  not  necessarily  restore  a  Sverdrup  type  balance 
in  the  interior  of  the  basin.  Figure  2.16  shows  a  sequence  of  time-averaged  stream  func¬ 
tion  and  potential  vorticity  fields  for  <5/  between  0.055  and  0.1725.  The  duration  of  the 
averaging  period  ranged  from  2500  nondimensional  time  units  to  10000  nondimensional 
units.  This  corresponds  to  averaging  periods  of  4  to  16  years  using  the  dimensional  scales 
given  in  Section  2.2.  As  the  forcing  increases  (increasing  6[)  the  western  intensification  of 
the  solution  decreases  despite  the  fact  that  the  solutions  do  not  become  highly  inertial. 
Only  close  to  the  marginal  stability  curve  does  the  time  averaged  solution  retain  some 
degree  of  western  intensification. 
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Figure  2.16:  Time-averaged  stream-function  (top  row),  and  potential  vorticity  field  (hot 
tom  row),  for  an  increasing  sequence  of  Sr-  Other  parameters  are  held  fixed  at  5q  —  0.0 


Chapter  3 


Reduced  Gravity  Model 

3.1  Model  Formulation 

In  this  chapter  we  will  deal  with  the  reduced  gravity  quasigeostrophic  vorticity  equation 
which  can  be  written  in  terms  of  the  interface  anomaly,  h,  as  follows 

(V^  -  +  (3h^  =  V^/i)  -  rV^h  -  +  ^^curlf.  (3.1) 

/o  Pog'H 

The  model  describes  the  time  evolution  of  the  interface  anomaly  between  two  immiscible, 
homogeneous  layers  of  fluid  of  slightly  different  densities.  The  lower  layer  is  assumed  to 
be  infinitely  deep  and  at  rest.  The  mathematical  model  is  the  same  as  that  described  in 
McCalpin  (1996)  and  McCalpin  and  Haidvogel  (1996). 

The  discretization  of  the  model  is  achieved  via  second-order  finite  difference  approx¬ 
imations,  with  a  horizontal  grid-spacing  of  20  km.  The  time  derivative  is  approximated 
using  Scheme  A  described  on  page  153  of  Haltiner  and  Williams  (1980)  with  a  time  step 
of  0.5  days. 

The  curl  of  the  wind  stress  is  given  by 

OtT  7/  7/1 
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Figure  3.1:  Wind-stress  curl  profile.  The  dashed  curve  {a  =  0.0)  is  exactly  anti¬ 
symmetric.  The  solid  curve  (a  =  0.05)  is  non-symmetric  with  the  peak  of  the  wind-stress 
curl  in  the  northern  gyre  decreased  by  5%  and  the  peak  of  the  wind-stress  curl  in  the 
southern  gyre  increased  by  5%,  relative  to  the  anti-symmetric  profile. 


The  parameter  a.  controls  the  north-south  asymmetry  of  the  wind-stress  curl.  Fig¬ 
ure  3.1  shows  the  wind-stress  curl  profile  for  a  =  0.0  and  a  =  0.05.  The  domain  of 
integration  is  a  rectangular  ocean  basin  confined  to  a  region  given  by 


V  -  {{x,y)  I  0  km  <x  <  3600  km,  and  0  km  <y  <  2800  km.}.  (3.3) 

At  the  basin  walls,  the  following  boundary  conditions  are  applied 

h  —  c(t),  V^/i  =  0,  and  =  0.  (3.4) 


The  quantity  c{t)  is  chosen  such  that  the  total  mass  in  the  upper  layer  remains  constant. 
This  is  achieved  by  solving 


' 

V^ht,  —  =0,  in  X> 

hb  =  1,  on  dV 

and 

'  -  'fho  =  V^h  -  in  V 

< 

ho  =  0,  on  dV 


(3.5) 


(3.6) 


64 


and  then  choosing  c{t)  such  that 


j I  h{x,  y,  t)dxdy  =  JJ  ho{x,  y,  t)  +  c{t)hb{x,  y)  dx  dy  =  0.  (3.7) 

The  numerical  values  for  the  parameters  are  the  following: 


North-south  extent  of  basin 
East-west  extent  of  basin 
Upper  layer  thickness 
Latitude  of  southern  basin  wall 
Coriolis  parameter 
Differential  rotation  of  the  earth 
Standard  density 
Reduced  gravity 
Rossby  radius  of  deformation 
Reciprocal  of  R^ 

Strength  of  wind-stress  control  parameter 
Wind-stress  profile  control  parameter 
Inter-facial  damping  coefficient 
Biharmonic  viscosity  coefficient 


Ly  =  2800  km, 

Lj.  =  3600  km, 

H  —  600  m, 

00  =  7r/6, 

/o  =  1.4544  X  lO-'‘cos(0o)  s-^ 
/?  =  2.2829  X  cos(0o) 

Po  =  1.027  X  10^  kg  m-^ 
g'  =  0.02  m  s-^, 

Rrf  =  47.636  km, 

7^  =  (47636)“^  m“^. 

To  =  0.05  Nm~^, 
a  =  0.05, 
r  —  10“^  s“^, 

Ab  =  8  X  101°  m'l  s-i. 


(3.8) 


3.2  Phenomenology 


To  reproduce  the  phenomenology  described  by  McCalpin  and  Haidvogel  (1996),  the 
model  was  run  for  more  than  1800  years,  saving  the  interface  height  field  every  5  days. 
Figure  3.2(a),  shows  a  1500  year  time  series  of  the  basin  integrated  total  energy,  and 
Figure  3.2(b)  shows  the  corresponding  time  series  of  the  eddy  energy.  It  can  be  seen  that 
periods  of  high  eddy  activity  are  associated  with  periods  where  the  model  is  in  a  state 
of  low  total  energy  and  conversely  that  periods  of  low  eddy  activity  are  associated  with 
periods  of  high  total  energy. 


65 


In  computing  the  time-series,  the  total  energy  was  expressed  as 


TE[h{t)]  =  PE[h{t)]  KE[h{t)l 


which  is  the  sum  of  the  potential  energy, 


and  the  kinetic  energy. 


PE[h{t)]  =  ^Pog'  J f  h^{t)dxdy, 
KE[h{t)]  =  j f  dx  dy, 


(3.10) 


where 


(3.11) 


(3.12) 


are  the  geostrophic  velocities. 

The  time  series  for  the  eddy  energy  was  obtained  by  subtracting  from  the  time-series 
of  total  energy  the  time-series  of  the  energy  computed  from  the  field  which  was  first 
low-pass  filtered  using  a  r  =  6  months  running  average. 


_ _  rt+T /2 

h  (t)  =  I  hit)  dt. 

Jt+T/2 


(3.13) 


The  eddy  or  high-frequency  energy  is  defined  to  be 


EE[h{t)]  =  TE[h{t)]  -  TE[h 


(3.14) 


In  Figure  3.3(a)  we  show  the  power  density  spectrum  for  the  potential  energy.  We 
can  see  that  the  spectrum  is  red  up  to  periods  of  50  years.  The  bend  in  the  spectrum  at 
periods  of  7  to  8  months  is  associated  with  the  preferred  period  of  eddy  jet  interactions. 
Figure  3.3(b)  shows  the  potential  energy  power  density  spectrum  multiplied  by  the 
frequency.  We  can  see  that  most  of  the  energy  lies  in  a  band  between  periods  of  10  years 
and  100  years.  Figure  3.4(a)  shows  the  kinetic  energy  power  density  spectrum.  The 
spectrum  is  red  up  to  periods  of  50  years,  and  is  white  for  longer  period.  There  is  also 
a  broad  peak  between  periods  of  6  months  to  2  years.  These  features  are  difficult  to  see 
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Frequency  X  Power  Density 


X  10^^  Potential  Energy  Spectrum 


Figure  3.3:  (a)  Potential  energy  power  density  spectrum  with  95%  confidence  intervals 
(b)  Frequency  x  potential  energy  power  density  spectrum. 


Figure  3.4:  (a)  Kinetic  energy  power  density  spectrum  with  95%  confidence  intervals 
(b)  Frequency  x  kinetic  energy  power  density  spectrum. 
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because  of  the  scale  of  the  ordinate,  but  are  nonetheless  significant.  Figure  3.4(b)  The 
power  density  spectrum  for  the  kinetic  energy  multiplied  by  the  frequency  (Figure  3.4(b)) 
shows  that  most  of  the  kinetic  energy  variability  occurs  at  periods  between  3  months  and 
100  years,  and  shows  the  peaks  between  6  months  and  2  years,  between  2  years  and  8 
years.  For  periods  longer  than  10  years  the  power  density  decreases. 

The  model  exhibits  low-frequency  variability  associated  with  irregular  transitions  be¬ 
tween  diflferent  dynamical  regimes  that  persist  for  extended  lengths  of  time.  McCalpin 
and  Haidvogel  (1996)  describe  these  regimes  as  follows:  a  high  energy  state  with  a  large 
jet  penetration  scale  and  weak  eddy/ring  formation,  a  low  energy  state  with  a  weakly 
penetrating  jet  and  strong  eddy/ring  generation,  and  an  intermediate  or  medium  energy 
state  with  intermediate  jet  penetration  and  modest  eddy/ring  formation.  They  further 
describe  these  regimes  as  having  time-averaged  interface  height  anomaly  fields  and  eddy 
energy  fields  which  are  remarkably  similar  from  event  to  event,  while  being  clearly  distinct 
from  one  regime  to  the  other.  Figure  3.5  is  reproduced  from  McCalpin  and  Haidvogel 
(1996)  and  shows  contour  plots  of  the  averaged  interface  height  anomaly  for  a  high, 
medium  and  low  energy  regime. 

3.3  Multiple  Equilibria 

For  finite  dimensional  dissipative  dynamical  systems  all  solutions  converge  as  t  oo 
to  a  complicated  set  called  the  global  attractor.  It  is  sometimes  the  case  that  unstable 
fixed  points  of  the  dynamical  system  lie  close  to,  or  are  contained  in  the  global  attractor 
(as  in  the  case  of  the  Lorenz  attractor;  Lorenz  (1963)).  When  this  is  the  case,  there  is 
the  possibility  that  the  recurrent  regimes  might  be  identified  with  the  model’s  unstable 
fixed  points.  In  general  there  is  a  stable  and  an  unstable  manifold  associated  with  each 
fixed  point.  The  unstable  manifold  is  the  nonlinear  extension  of  the  solution’s  unstable 
tangent  space  (i.e.  the  space  spanned  by  the  unstable  eigenvectors  of  the  model  linearized 
about  the  fixed  point).  Almost  all  trajectories  which  begin  in  a  neighborhood  of  the  fixed 
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Figure  3.5:  Time  averaged  interface  anomaly  fields  from  McCalpin  and  Haidvogel  (1996). 
Upper  left  is  for  a  200  year  record.  Upper  right  is  for  a  medium  energy  period.  Lower 
left  is  for  a  high  energy  period  and  Lower  right  is  for  a  low  energy  period.  All  frames 
employ  the  same  contour  interval  of  20  m.  The  axis  are  horizontal  distances  in  km. 


point  are  eventually  expelled  along  the  unstable  manifold.  These  trajectories  can  however 
sometimes  return  to  the  neighborhood  of  the  fixed  point  by  following  a  trajectory  close 
to  the  stable  manifold.  In  this  sense  the  unstable  fixed  points  act  to  “steer”  the  time 
dependent  behavior  (Legras  and  Ghil  1985). 

To  test  the  possibility  that  unstable  fixed  points  do  indeed  act  to  steer  the  model 
trajectories  in  McCalpin  and  Haidvogel’s  simulations,  an  arclength  continuation  method 
was  used  to  search  for  the  fixed  points  of  the  system.  By  first  setting  a  to  zero  and 
gradually  decreasing  the  biharmonic  viscosity  coefficient,  Ai,,  the  anti-symmetric  solution 
branch  was  mapped  from  a  linear  viscous  regime  to  the  non-linear  regime  of  McCalpin 
and  Haidvogel’s  simulation.  Recall  that  a  is  the  parameter  controlling  the  asymmetry  of 
the  wind-stress  curl.  Initially  the  wind-stress  curl  was  kept  anti-symmetric  by  setting  a 
to  zero  so  that  any  pitchfork  bifurcation  structure  leading  to  multiple  equilibria  will  not 
be  destroyed  (Jiang  et  al.  (1995),  Cessi  and  lerley  (1995)).  For  a  different  from  zero  the 
branches  which  would  be  connected  at  the  pitchfork  bifurcation  point  for  a  =  0,  become 
disconnected  and  cannot  all  be  found  by  continuously  varying  the  viscosity  coefficient. 
Once  a  solution  on  each  of  the  distinct  branches  was  found,  the  branches  were  extended  to 
values  of  a  different  from  zero.  Furthermore,  to  reduce  the  computational  costs,  the  fixed 
points  were  first  calculated  on  a  coarse  grid.  The  need  to  resolve  the  sharp  gradients  in 
the  potential  vorticity  field  along  the  western  wall,  and  in  the  inter-gyre  regions  requires 
the  use  of  a  stretched  grid.  Figure  B.l  in  Appendix  B  shows  the  coarse  computational 
grid  with  grid  points  concentrated  along  the  western  wall  and  the  basin  center.  The  grid 
spacing  varies  in  the  X  direction  from  a  minimum  of  30  km  along  the  western  wall  to 
a  maximum  of  90  km  along  the  eastern  wall,  and  in  the  Y  direction,  the  grid  spacing 
varies  from  a  minimum  of  30  km  in  the  center  of  the  basin  to  a  maximum  of  90  km  along 
the  northern  and  southern  walls.  The  resolution  used  to  compute  each  solution  was  then 
gradually  increased  to  the  final  value  of  Ax  =  Ay  =  20  km,  -  the  same  resolution  used  in 
the  time-dependent  simulation.  Also  see  Section  3.3.12  for  a  discussion  of  the  difficulties 
associated  with  interpolating  steady  state  solutions  onto  successively  finer  grids. 
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3.3.1  Anti-symmetric  Solutions:  (a  =  0.0) 

Starting  from  the  linear  viscous  Munk  like  solution,  the  biharmonic  viscosity  coefficient, 
A/,,  was  gradually  reduced  to  allow  the  steady  state  solution  to  become  progressively  more 
nonlinear.  Figure  3.6  shows  the  maximum  transport  in  the  jet  for  the  anti-symmetric 
solution  as  a  function  of  the  viscosity  parameter.  It  also  shows  the  location  of  the  bi¬ 
furcation  points  about  which  will  be  further  explained  below.  Recirculation  cells  appear 
when  the  biharmonic  viscosity  coefficient,  Ab,  is  reduced  to  1.5  x  10^^  m^s“^.  As  it  is  fur¬ 
ther  reduced,  the  recirculation  cells  intensify  and  expand  eastward.  With  the  increasing 
intensity  of  the  recirculation  cells,  there  is  a  sharp  increase  in  the  maximum  transport 
across  the  jet  until  a  maximum  of  42  Sv  is  reached  for  Ab  =  1.3  x  10^^  m^s“F  When  this 
maximum  in  the  transport  is  reached,  the  recirculation  cells  extend  600  km  eastward  into 
the  basin  interior.  A  subsequent  reduction  in  the  viscosity  causes  the  intensity  of  the 
recirculation  cells  io  decrease,  but  their  eastward  extent  continues  to  increase  in  a  mono¬ 
tonic  fashion  until  a  saddle  node  bifurcation  point  is  reached  at  Ab  =  6.0  x  10^  m‘*s“F 
The  low  nose  point  associated  with  this  bifurcation  is  labeled  NPL  in  Figure  3.6.  At 
this  point  the  recirculation  cells  extend  up  to  180  km  west  of  the  eastern  wall.  To  contin¬ 
uously  follow  the  anti-symmetric  solution,  the  viscosity  must  be  increased  from  the  low 
nose  point,  NPL,  up  to  the  high  nose  point,  NPH,  at  Ab  =  7.6  x  10^  m'*s“^  Once  this 
second  saddle  node  bifurcation  is  reached,  the  recirculation  cells  extend  fully  across  the 
basin  from  the  western  wall  to  the  eastern  wall.  A  subsequent  decrease  in  the  viscosity 
coefficient  causes  the  recirculation  cells  to  expand  in  size  in  the  north-south  direction,  and 
to  intensify.  This  intensifying  of  the  recirculation  cells  causes  the  maximum  transport  to 
increase  again. 

In  addition  to  the  two  saddle  node  bifurcations  described  above,  eight  successive 
symmetry  breaking  pitchfork  bifurcations  occur  as  the  viscosity  is  decreased.  They  are 
labeled  PFa,  PFg,  PFc,  PFu,  PFe,  PFp,  PFq,  and  PFh  in  Figure  3.6.  In  Figures  3.7 
to  3.16  the  steady  state  solution  at  each  of  the  ten  bifurcation  points  (8  pitchforks  -I-  2 
saddle  nodes),  as  well  as  the  corresponding  null  eigen-mode  (i.e.  the  one  with  the  zero 
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Maximum  Transport 


Figure  3.6:  Plot  of  the  maximum  transport  across  the  inter-gyre  jet  for  the  anti-symmetric 
solutions  as  a  function  of  the  biharmonic  viscosity  coefficient,  Ai,.  Also  shown  is  the 
location  of  the  pitchfork  bifurcation  points,  PFa,  PFb,  PFc,  PFd,  PFe,  PFp,  PFq, 
and  PFf{,  as  well  as  the  location  of  two  saddle  node  bifurcation  points  NPL  and  NPH. 


eigenvalue)  are  contoured  in  the  order  in  which  they  occurred  as  the  pseudo-arclength 
was  increased  (see  Appendix  B  for  a  definition  of  pseudo-arclength).  At  bifurcation 
points,  one  of  the  eigenmodes  of  the  system  linearized  about  the  fixed  point  has  a  zero 
eigenvalue.  This  null  mode  captures  the  essential  difference  between  the  anti-symmetric 
fixed  point  at  the  bifurcation  point  and  the  new  equilibria  that  come  into  existence  as 
a  result  of  the  bifurcation.  In  the  following  sections  each  of  the  bifurcations  are  now 
discussed  in  turn.  The  following  sections  will  show  how  each  pitchfork  bifurcations  adds 
an  additional  half  meander  in  the  jet  separating  the  recirculations  cells.  The  meanders 
that  occur  up-stream  where  the  jet  is  strongest  have  the  longest  wave-length,  and  those 
that  occur  down-stream  have  shorter  wave-lengths.  This  suggests  that  the  meanders 
might  be  understood  in  terms  of  damped  stationary  Rossby  vaves.  Further  investigation 
along  these  lines  is  left  for  future  work. 

3.3.2  First  Bifurcation:  PFa 

The  first  pitchfork  bifurcation  point,  labeled  PFa,  {At,  =  4.1  x  10^^  m^s“^),  occurs  shortly 
after  the  recirculation  gyres  form,  when  their  eastward  extent  is  only  400  km  into  the  basin 
interior  (Figure  3.7).  The  associated  null  eigen-mode  is  symmetric  about  the  line  of  zero 
wind-stress  curl  (which  must  be  the  case  for  a  symmetry  breaking  pitchfork  bifurcation). 
The  mode’s  structure  consists  of  three  counter  rotating  gyres.  The  strongest  is  situated 
over  the  jet  axis  and  extends  550  km  into  the  basin  interior.  It  is  roughly  pear  shaped 
with  its  widest  extent  overlying  the  eastern  most  part  of  the  recirculation  gyres.  To  the 
north  and  south,  and  slightly  to  the  west,  it  is  flanked  by  counter  rotating  gyres.  When 
the  mode  is  superimposed  on  the  basic  state,  it  causes  one  of  the  recirculation  cells  to 
intensify  and  the  other  to  weaken  depending  on  the  sign  of  the  mode.  The  superposition 
also  causes  the  recirculation  cell  pattern  to  shift  either  north  or  south  with  the  weaker 
recirculation  cell  being  pulled  across  the  line  of  zero  wind-stress  curl  by  its  stronger 
counter  part. 
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Figure  3.7:  (a)  Contour  plot  of  the  interface  height  anomaly,  h,  (C.I.  =  20  m)  at  the 
pitchfork  bifurcation  point  PFa-  (b)  Contour  plot  of  the  null  mode  with  the  zero  eigen¬ 
value  responsible  for  the  bifurcation  (amplitude  and  sign  are  arbitrary),  (c)  Zoom  in  of 
interface  height  anomaly,  (d)  Zoom  in  of  the  null  mode. 
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3.3.3  Second  Bifurcation:  PFb 

The  second  pitchfork  bifurcation  point,  labeled  PFb,  {Ab  =  5.7  x  10^^  occurs 

after  the  recirculation  gyres  have  reached  their  maximum  intensity  and  have  begun  to 
weaken.  At  the  bifurcation  point,  the  recirculation  gyres  extend  eastward  slightly  more 
than  800  km  (Figure  3.8).  The  associated  null  eigen-mode  is  again  symmetric,  and  its 
structure  consists  of  three  counter  rotating  gyres.  It  has  a  maximum  situated  over  the 
jet  axis,  in  the  same  way  that  the  null  mode  did  for  the  bifurcation  point  PFa,  but  this 
time,  the  strongest  of  the  three  cells  which  is  situated  over  the  jet  axis  does  not  extend 
east  of  the  recirculation  gyres.  Its  eastern  most  edge  coincides  with  the  eastern  most  edge 
of  the  recirculation  gyres.  The  structure  of  the  mode  in  the  region  of  the  jet  axis  is  tear 
shaped,  and  the  zonal  position  of  its  maximum  is  nearly  coincident  with  the  maximum 
for  the  recirculation  gyres.  To  the  north  and  south,  it  is  flanked  by  two  weaker  counter 
rotating  gyres  which  wrap  around  the  eastern  edge  of  the  tear  shaped  structure.  When 
it  is  added  to  the  basic  state,  it  causes  one  recirculation  cell  to  intensify  and  the  other  to 
weaken,  depending  on  the  sign  of  the  mode.  It  also  causes  the  jet  to  turn  in  the  direction 
of  the  stronger  recirculation  cell  after  it  separates  from  the  western  wall  at  the  line  of 
zero  wind-stress  curl.  But  unlike  the  null  mode  at  PFA  which  caused  the  jet  to  turn 
fully  on  itself,  this  null  mode  causes  the  jet  to  turn  eastward  after  moving  away  from  line 
of  zero  wind-stress  curl,  and  only  then  does  it  fan  out  to  rejoin  the  Sverdrup  interior. 
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Figure  3.8:  (a)  Contour  plot  of  the  interface  height  anomaly,  h,  (C.I.  =  20  m)  at  the 
pitchfork  bifurcation  point  PFb-  (b)  Contour  plot  of  the  null  mode  with  the  zero  eigen¬ 
value  responsible  for  the  bifurcation  (amplitude  and  sign  are  arbitrary),  (c)  Zoom  in  of 
interface  height  anomaly,  (d)  Zoom  in  of  the  null  mode. 
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3.3.4  Third  Bifurcation:  PFc 

The  third  pitchfork  bifurcation  point,  labeled  PFc,  {At,  =  2.4  x  10^^  occurs  when 

the  recirculation  gyres  extend  1250  km  eastward  into  the  basin  interior  (Figure  3.9).  The 
associated  null  eigen-mode  again  consistis  of  three  counter  rotating  cells.  The  cell  that  is 
situated  over  the  jet  axis  consists  of  an  elongated  structure  with  two  local  maximums  at 
650  km  and  1250  km.  The  saddle  point  separating  the  maximums  is  situated  at  1100  km 
on  the  jet  axis,  so  that  the  “bump”  associated  with  the  first  maximum  is  1100  km  long, 
and  the  one  associated  with  the  second  maximum  is  only  400  km  long  since  this  mode  does 
not  extend  past  1500  km.  Like  the  null  mode  associated  with  the  bifurcation  PFa,  the 
null  mode  at  PFc  has  a  pear  shaped  structure  that  extends  eastward  of  the  region  where 
the  jet  begins  to  fan  out  to  rejoin  the  Sverdrup  interior.  The  other  two  cells  are  situated 
to  the  north  and  south  and  are  positioned  over  the  westward  recirculating  currents.  When 
the  mode  is  superimposed  on  the  basic  state  it  causes  one  of  the  recirculation  cells  to 
intensify  and  the  other  to  weaken  depending  on  the  sign  of  the  mode.  It  also  shifts  the 
jet  axis  either  north  or  south  in  the  direction  of  the  stronger  recirculation  gyre.  The 
mode  also  causes  a  meander  in  the  jet.  After  separating  from  the  coast  at  the  line  of  zero 
wind-stress  curl,  the  perturbed  jet  turns  in  the  direction  of  the  stronger  recirculation, 
then  turns  eastward  again,  and  finally  turns  back  to  the  line  of  zero  wind-stress  curl, 
before  fanning  out  to  rejoin  the  Sverdrup  interior. 
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Steady  state  at  PF^  =  2.4e+1 1 ) 
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Figure  3.9;  (a)  Contour  plot  of  the  interface  height  anomaly,  h,  (C.I.  =  20  m)  at  the 
pitchfork  bifurcation  point  PFc-  (b)  Contour  plot  of  the  null  mode  with  the  zero  eigen¬ 
value  responsible  for  the  bifurcation  (amplitude  and  sign  are  arbitrary),  (c)  Zoom  in  of 
interface  height  anomaly,  (d)  Zoom  in  of  the  null  mode. 
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3.3.5  Fourth  Bifurcation:  PFd 

The  fourth  pitchfork  bifurcation  point,  labeled  PFd,  {At,  =  1.2  x  10^^  occurs 

when  the  recirculation  gyres  extend  1750  km  eastward  into  the  basin  interior  (Fig¬ 
ure  3.10).  For  the  most  part,  the  associated  null  mode  is  similar  to  the  null  mode 
associated  with  the  bifurcation  point  PFc-  It  is  anti-symmetric  and  composed  of  three 
cells,  with  the  cell  overlying  the  jet  axis  having  two  relative  maxima.  Like  the  null  mode 
for  PFc,  the  first  maximum  is  slightly  to  the  east  of  the  center  of  a  1100  km  long  “bump” . 
The  eastern  part  of  the  mode  however  is  most  like  the  null  mode  associated  with  the  bi¬ 
furcation  PFb-  It  has  a  structure  which  is  tear  shaped  and  is  nearly  coincident  with 
the  location  of  maximum  transport  in  the  jet.  The  second  maximum  is  72  %  as  intense 
as  the  first.  Like  the  null  mode  for  PFb,  the  two  counter  rotating  gyres  wrap  around 
the  anomaly  that  is  situated  over  the  jet  axis.  Also  like  the  null  mode  associated  with 
PFb,  but  unlike  the  modes  associated  with  PF^  and  PFc,  this  mode’s  eastern  most 
edge  coincides  with  the  eastern  most  edge  of  the  basic  states’  recirculation  gyres.  When 
this  mode  is  superimposed  on  the  basic  state  it  causes  the  jet  to  shift  either  southward 
or  northward  depending  on  the  sign  of  the  mode.  It  also  causes  the  jet  to  meander  in  a 
similar  way  that  the  null  mode  of  PFc  forces  the  jet  to  meander,  except  that  the  jet  turns 
one  additional  time  away  from  the  line  of  zero  wind-stress  curl  before  turning  eastward 
and  fanning  out  to  rejoin  the  Sverdrup  interior. 
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Steady  state  at  PF^  (Aj^  =  1 .2e+1 1 ) 


Null  eigenmode  at  PF^ 


Steady  state  at  PF^  =  1 .2e+1 1 )  Null  eigenmode  at  PF^ 


Figure  3.10;  (a)  Contour  plot  of  the  interface  height  anomaly,  h,  (C.I.  =  20  m)  at  the 
pitchfork  bifurcation  point  PFo-  (b)  Contour  plot  of  the  null  mode  with  the  zero  eigen¬ 
value  responsible  for  the  bifurcation  (amplitude  and  sign  are  arbitrary),  (c)  Zoom  in  of 
interface  height  anomaly,  (d)  Zoom  in  of  the  null  mode. 
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3.3.6  Fifth  Bifurcation:  PFe 


The  fifth  pitchfork  bifurcation  point,  labeled  PFe,  (A  =  6  x  10^“  rn'^s"’^),  occurs  when 
the  recirculation  gyres  extend  2100  km  eastward  into  the  basin  interior  (Figure  3.11).  The 
null  mode  associated  with  this  bifurcation  point  has  a  dominant  structure  with  three  local 
maximums  overlying  the  jet  axis.  The  saddle  points  separating  the  three  “bumps”  are 
situated  at  1150  km  and  2100  km.  The  third  and  eastern  most  bump  extends  from  the 
second  saddle  point  up  to  2400  km.  The  first  “bump”  is  1150  km  long.  The  second 
is  1050  km  long,  and  the  third  is  300  km  long.  The  second  maximum  has  an  intensity 
which  is  82%  of  the  first,  and  the  third  has  an  intensity  which  is  30%  of  the  first.  To  the 
north  and  south,  the  null  mode  has  two  counter  rotating  gyres  which  extend  from  the 
western  wall  up  to  the  same  zonal  position  as  that  of  the  maximum  for  the  third  bump. 
These  cells  have  maximums  near  2100  km,  the  same  zonal  position  as  that  of  the  second 
saddle  point  separating  the  second  and  third  bump.  The  eastern  structure  of  this  mode 
is  similar  to  the  null  modes  associated  with  the  bifurcation  points  PF^  and  PFc-  Like 
the  modes  for  those  previous  bifurcations,  the  eastern  most  structure  is  pear  shaped  and 
extends  to  the  east  of  the  recirculation  gyres.  When  this  mode  is  added  to  the  basic  state 
its  effect  is  similar  to  that  of  the  null  mode  associated  with  PFc  except  that  the  jet  has 
an  additional  meander. 
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Figure  3.11:  (a)  Contour  plot  of  the  interface  height  anomaly,  h,  (C.I.  :=  20  m)  at  the 
pitchfork  bifurcation  point  P Fe-  (b)  Contour  plot  of  the  null  mode  with  the  zero  eigen¬ 
value  responsible  for  the  bifurcation  (amplitude  and  sign  are  arbitrary),  (c)  Zoom  in  of 
interface  height  anomaly,  (d)  Zoom  in  of  the  null  mode. 


3.3.7  Sixth  Bifurcation:  PFp 

The  sixth  pitchfork  bifurcation  point,  labeled  PFpX^b  =  2.8  x  10^*^  m^s“^) ,  occurs  when 
the  recirculation  gyres  extend  approximately  2500  km  eastward  into  the  basin  interior 
(Figure  3.12).  The  null  mode  associated  with  this  bifurcation  has  a  similar  structure 
to  the  mode  associated  with  PFq,  except  that  three  instead  of  two  local  maxima  exist 
on  the  cell  which  is  situated  over  the  jet  axis.  The  first  saddle  point  separating  the 
maxima  is  situated  at  1200  km  and  the  second  at  2150  km.  The  eastern  structure  of 
the  mode  has  the  familiar  tear  shaped  structure  also  present  for  the  null  modes  of  PFb 
and  PFd-  Like  those  previous  bifurcations  the  tear  shaped  structure  does  not  extend 
past  the  recirculation  gyres  of  the  basic  state.  The  mode  also  has  the  familiar  counter 
rotating  cells  to  the  north  and  south  which  wrap  around  the  most  eastern  extent  of  the 
structure. 
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steady  state  at  PF^  (A^  =  2.8e+10)  Null  eigenmode  at  PFp 
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Figure  3.12:  (a)  Contour  plot  of  the  interface  height  anomaly,  h,  (C.I.  =  20  m)  at  the 
pitchfork  bifurcation  point  P Pp.  (b)  Contour  plot  of  the  null  mode  with  the  zero  eigen¬ 
value  responsible  for  the  bifurcation  (amplitude  and  sign  are  arbitrary),  (c)  Zoom  in  of 
interface  height  anomaly,  (d)  Zoom  in  of  the  null  mode. 
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3.3.8  Seventh  Bifurcation:  PFq 

The  seventh  pitchfork  bifurcation  point,  labeled  PFg,  (A  =  1.3  x  10^“  occurs 

when  the  recirculation  gyres  extend  approximately  3000  km  eastward  into  the  basin  inte¬ 
rior.  The  null  mode  associated  with  this  bifurcation  is  similar  to  the  null  mode  associated 
with  the  bifurcation  PFe,  except  that  this  mode  causes  one  additional  meander  in  the 
jet  when  it  is  added  to  the  basic  state.  Also  different  from  all  the  previous  modes,  the 
bumps  on  the  cell  overlying  the  jet  axis  are  separated  by  small  valleys  instead  of  saddle 
points.  A  further  difference  is  that  this  mode  causes  the  meanders  in  the  jet  to  overshoot 
the  line  of  zero  wind-stress  curl  when  it  is  added  to  the  basic  state.  This  difference  might 
be  due  to  the  fact  that  the  null  mode  feels  the  presence  of  the  eastern  wall  which  lies 
less  than  100  km  to  the  eastern  edge  of  the  last  cell  on  the  jet  axis.  The  intensity  of 
the  co-rotating  cells  on  the  jet  axis  decreases  downstream  such  that  the  intensity  of  the 
second  cell  is  72%  that  of  the  first,  the  intensity  of  the  third  is  68%  of  the  first,  and 
that  of  the  fourth  is  32%  the  intensity  of  the  first.  The  intensity  of  the  weaker  counter 
rotating  cells  increases  downstream  such  that  the  intensity  of  the  first,  second,  third  and 
fourth  “valley”  is  7%,  8%,  12%  and  15%  the  intensity  of  the  strongest  “bump”. 
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steady  state  at  PF^  =  1 .3e+1 0)  Null  eigenmode  at  PF^ 


Figure  3.13:  (a)  Contour  plot  of  the  interface  height  anomaly,  h,  (C.I.  =  20  m)  at  the 
pitchfork  bifurcation  point  P Fq.  (b)  Contour  plot  of  the  null  mode  with  the  zero  eigen¬ 
value  responsible  for  the  bifurcation  (amplitude  and  sign  are  arbitrary),  (c)  Zoom  in  of 
interface  height  anomaly,  (d)  Zoom  in  of  the  null  mode. 
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3.3.9  Eight  and  Ninth  Bifurcations:  NPL  and  NPH 

The  next  two  bifurcation  points  are  saddle  nodes  associated  with  the  interaction  of  the 
jet  and  the  eastern  wall.  The  first,  NPL,  is  the  low  nose  point  which  occurs  when  the 
jet  reaches  the  eastern  wall  at  Ab  =  6.0  x  10^  m^s"^  (Figure  3.14).  The  associated  null 
eigen-mode  is  anti-symmetric  like  the  basic  state,  and  consists  of  two  counter  rotating 
cells  trapped  near  the  eastern  wall  in  the  region  where  the  streamlines  in  the  jet  fan 
out.  The  second  saddle  node  bifurcation  point  (Figure  3.15)  is  labeled  NPH,  and  is 
the  high  nose  point.  It  occurs  for  Ab  =  7.6  x  10^  m^s“^.  The  associated  eigen-mode  is 
anti-symmetric  and  consists  of  two  counter  rotating  cells  which  are  most  intense  near 
the  eastern  wall,  but  which  extend  westward  up  to  the  western  wall.  When  the  mode  is 
added  to  the  basic  state  it  causes  the  recirculation  gyres  to  expand  or  recede  northward 
and  southward,  depending  on  the  sign  of  the  mode. 
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steady  state  at  NPL  (A^  =  6.0e+09) 
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Figure  3.14:  (a)  Contour  plot  of  the  interface  height  anomaly,  h,  (C.I.  =  20  m)  at  the  low 
nose  point  NPL.  (b)  Contour  plot  of  the  null  mode  with  the  zero  eigen-value  responsible 
for  the  bifurcation  (amplitude  and  sign  are  arbitrary),  (c)  Zoom  in  of  interface  height 
anomaly,  (d)  Zoom  in  of  the  null  mode. 
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Figure  3.15:  (a)  Contour  plot  of  the  interface  height  anomaly,  h,  (C.I.  =  20  m)  at  the  high 
nose  point  NPH.  (b)  Contour  plot  of  the  null  mode  with  the  zero  eigen-value  responsible 
for  the  bifurcation  (amplitude  and  sign  are  arbitrary),  (c)  Zoom  in  of  interface  height 
anomaly,  (d)  Zoom  in  of  the  null  mode. 


3.3.10  Tenth  Bifurcation:  PFh 


The  eighth  pitchfork  bifurcation  point,  labeled  PFh  ,  (A*  =  1.3  x  10^”  occurs 

after  the  jet  extends  across  the  basin,  and  the  recirculation  cells  have  begun  to  expand 
northward  and  southward  (Figure  3.16).  The  null  mode  associated  with  the  bifurcation 
consists  of  4  co-rotating  cells  separated  by  very  weak  counter  rotating  cells  which  all  lie 
on  the  jet  axis.  The  three  valleys  separating  the  four  local  maximums  are  situated  near 
1300  km,  2300  km,  and  3200  km.  The  second  maximum  is  60%  as  strong  as  the  first, 
the  third  is  50%  as  strong  as  the  first  and  the  third  is  35%  as  strong  as  the  first.  When 
the  mode  is  superimposed  on  the  basic  state  it  causes  the  jet  to  shift  either  northward 
or  southward  depending  on  the  sign  of  the  mode.  It  also  causes  the  jet  to  have  four 
meanders  of  length  1300  km,  1000  km,  900  km  and  400  km. 
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Figure  3.16:  (a)  Contour  plot  of  the  interface  height  anomaly,  h,  (C.I.  =  20  m)  at  the 
pitchfork  bifurcation  point  PFh-  (b)  Contour  plot  of  the  null  mode  with  the  zero  eigen¬ 
value  responsible  for  the  bifurcation  (amplitude  and  sign  are  arbitrary),  (c)  Zoom  in  of 
interface  height  anomaly,  (d)  Zoom  in  of  the  null  mode. 


3.3.11  Non-symmetric  Solutions:  (o  =  0.0) 

Four  of  the  eight  symmetry  breaking  pitchfork  bifurcations  discussed  above  occurred  for 
values  of  Ab  less  than  8  x  10^°m^s"\  so  that  for  Ab  =  8  x  10^“m'‘s-^  (the  value  used 
by  McCalpin  and  Haidvogel  (1996))  and  a  =  0.0,  there  exists  9  distinct  steady  state 
solutions  -  one  which  is  anti-symmetric  and  four  pairs  of  solutions  which  have  the  same 
symmetry  as  the  model,  ij){x,y  -  Ly/2)  =  -tp{x,  -{y  -  Ly/2)). 

Figure  3.17  shows  a  plot  of  energy  as  a  function  of  the  biharmonic  viscosity  coef¬ 
ficient  Ab.  The  plot  also  shows  the  non-symmetric  branches  which  have  bifurcated  at 
the  symmetry  breaking  pitchfork  bifurcations  PFa,  PFb,  PFc,  and  PFd-  Since  the 
wind-stress  curl  profile  for  a  =  0  is  exactly  anti-symmetric,  the  members  of  each  pair  of 
non-symmetric  equilibria  have  the  same  energy  and  thus  fall  on  overlapping  curves  which 
are  labeled  i^A^A' ,  B&cB',  C^C ,  and  D&cD').  The  range  between  Ab  =  lO^^m^/s  and 
Ab  =  lO^^m^/s  in  which  the  energy  of  the  flow  remains  essentially  constant  corresponds 
to  the  range  of  parameters  where  the  flow  is  essentially  linear  with  Sverdrup  balance 
everywhere  in  the  interior  of  the  basin  except  for  a  thin  western  boundary  layer  of  thick¬ 
ness  Sff  =  (Aj,/^)^/^.  For  Ab  near  lO^^m^/s,  recirculation  cells  form  in  the  region  where 
the  western  boundary  currents  from  the  sub-tropical  and  sub-polar  gyres  meet.  As  the 
biharmonic  diffusivity  is  further  decreased,  the  total  energy  of  the  anti-symmetric  branch 
increases  rapidly.  This  rapid  increase  in  energy  is  accompanied  by  a  rapid  increase  in 
the  zonal  extent  of  the  recirculation  cells.  As  the  recirculation  cells  continuously  expand 
in  the  zonal  direction  they  allow  the  successive  pitchfork  bifurcations  to  occur.  Each 
pair  of  new  equilibria  has  an  additional  half  meander  in  the  part  of  the  jet  separating 
the  counter  rotating  recirculation  gyres.  Thus  equilibria  A  and  A'  simply  turn  back  on 
themselves,  one  to  the  north  and  the  other  to  the  south.  Equilibria,  B  and  B'  first  turn 
north  or  south,  and  then  turn  seaward  before  fanning  out  to  rejoin  the  Sverdrup  interior. 
Equilibria  C  and  C  add  half  an  additional  meander,  and  so  on  for  D  and  D' .  In  contrast 
to  the  anti-symmetric  branch,  the  energy  level  for  the  non-symmetric  branches  A&A', 
BScB',  C'&C",  and  Dk.D'  decreases  or  remains  nearly  constant  (Figure  3.17).  For  these 
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Figure  3.17:  Plot  of  total  energy  divided  by  (TE)  for  a  =  0  as  a  function  of  the  lateral 
diffusion  parameter  Ab.  The  dashed  vertical  line  at  Ab  =  8*®  gives  the  value  of  lateral 
diffusion  used  in  the  time  dependent  simulation.  The  circles  indicate  the  bifurcation 
points  for  the  anti-symmetric  branch.  The  branches  labeled  AhA' ,  BSzB',  CfeC",  and 
Dk.D'  are  for  non-symmetric  solutions.  The  anti-symmetric  branch  is  the  one  that  exists 
for  all  parameter  values.  The  anti-symmetric  branch  would  become  £/'  at  =  8  x  10^° 
if  a  =  0.05. 
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non-symmetric  solutions,  large  inter-gyre  fluxes  of  vorticity  allow  the  solutions  to  remain 
balanced.  The  recirculation  cells  do  not  need  to  expand  in  size  to  allow  the  explicit  dissi¬ 
pation  to  remove  the  excess  vorticity  which  would  otherwise  have  been  dissipated  in  the 
western  boundary  layer,  or  more  efficiently  across  the  jet  axis  by  a  larger  lateral  diffusion 
parameter.  The  gyre  integrated  vorticity  budget  is  further  discussed  in  Section  3.3.13 
below. 

3.3.12  Non-symmetric  Solutions:  (a  =  0.05) 

As  a  was  increased  from  0.0  to  0.05  the  branch  labeled  A  in  Figure  3.17  turned  (i.e. 
underwent  a  saddle  node  bifurcation)  before  a  could  be  increased  to  0.05.  All  the  other 
branches,  however,  could  be  traced  continuously  from  a  =  0.0  to  0.05  .  Thus,  8  distinct 
steady  state  solutions  were  found  for  the  parameter  values  used  in  the  time-dependent 
simulation.  Contour  plots  of  the  layer  thickness  and  potential  vorticity  field  for  each 
of  the  equilibrium  states  computed  on  a  uniform  grid  with  20  km  resolution  are  shown 
in  Figures  3.18  through  3.25.  Most  of  these  solutions  were  obtained  on  the  fine  grid 
by  successively  interpolating  the  coarse  resolution  solution  onto  a  grid  with  a  few  more 
grid  points,  performing  several  Newton  correction  steps,  and  repeating  the  process  until 
the  desired  resolution  was  obtained.  However,  for  the  non-symmetric  solutions  A,  A', 
B,  and  B',  which  are  on  branches  which  have  turning  points  for  values  of  Ab  less  than 
8  X  10^°  m^s  \  (see  Figure  3.17),  Newton’s  method  did  not  always  converge  even  when 
only  one  grid  point  was  added  at  a  time.  For  those  branches  solutions  at  higher  values  of 
the  biharmonic  viscosity  had  to  be  interpolated  onto  successively  finer  grids  for  Newton’s 
method  to  converge.  The  high  resolution  branches  could  then  be  traced  back  to  the  lower 
values  of  the  viscosity. 

The  solution  labeled  E'  is  from  the  solution  branch  which  would  be  anti-symmetric 
if  the  wind-stress-curl  profile  was  exactly  anti-symmetric,  i.e.  if  a  =  0.0.  We  label  it  E' 
because  for  a  5^  0,  the  pitchfork  bifurcation  is  destroyed  and  the  branch  connects  to  the 
E'-branch  without  going  through  a  bifurcation  point.  The  pair  of  equilibria  labeled  D 
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Figure  3.18:  Potential  vorticity  field  (top,  C.I.  3.1  x  10~®  )  and  interface  height  anomaly 
(bottom,  C.I.  20  m)  for  fixed  point  A. 
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Figure  3.19:  Potential  vorticity  field  (top,  C.I.  3.1  x  10'^  )  and  interface  height  anomaly 
(bottom,  C.I.  20  m)  for  fixed  point  B. 
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Figure  3.20:  Potential  vorticity  field  (top,  C.I.  3.1  x  10  ^  )  and  interface  height  anomaly 
(bottom,  C.I.  20  m)  for  fixed  point  B'. 
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C:  Potential  Vorticity  (q) 


Figure  3.21:  Potential  vorticity  field  (top,  C.I.  3.1  x  10"®  )  and  interface  height  anomaly 
(bottom,  C.I.  20  m)  for  fixed  point  C. 


and  D'  are  from  the  branches  which  bifurcated  at  a  value  of  the  viscosity  closest  to  that 
used  in  the  time-dependent  simulation.  The  equilibria  C  and  C  bifurcated  at  a  larger 
value  of  viscosity,  and  the  pair  B  and  B'  bifurcated  at  a  yet  larger  value  of  viscosity. 
Finally  equilibrium  A  is  from  one  of  the  two  branches  that  bifurcated  at  the  largest  value 
of  viscosity.  The  equilibrium  A'  is  the  one  that  could  not  be  traced  continuously  to  a 
value  of  a  =  0.05. 

The  major  difference  in  the  structure  of  the  flow  fields  among  the  various  pairs  of 
equilibria  lies  in  the  extent  of  penetration  of  the  jet  into  the  interior,  and  in  the  degree 
to  which  the  streamlines  in  the  jet  meander  before  rejoining  the  Sverdrup  interior.  Both 
these  features  will  be  revisited  in  Section  3.3.15  where  we  examine  the  q  —  '^  relationship 
in  the  jet.  At  the  bifurcation  points,  all  the  flow  fields  are  exactly  anti-symmetric.  The 
amount  of  meandering  and  wave  activity  increases  continuously  for  the  non-symmetric 
branches  as  the  viscosity  is  continuously  decreased  from  its  value  at  the  bifurcation 
point.  The  penetration  scale  of  the  jet  on  the  other  hand  increases  only  slightly  for 
the  non-symmetric  branches  as  the  viscosity  coefficient  is  similarly  decreased.  Compare 
Figure  3.7  which  show  the  jet  which  penetrates  only  400  km  at  the  bifurcation  point 
PFa,  and  Figure  3.18  which  shows  the  non-symmetric  equilibria,  A,  after  the  viscosity 
has  been  reduced  from  4.1  x  10^^  m^s“^  to  8.0  x  10^°  m^s“^  to  see  that  the  jet  penetrates 
no  more  than  450  km  before  it  begins  to  fan  out.  Similarly  equilibrium  B  penetrates 
800  km  at  its  bifurcation  point  {A^  —  5.7  x  10^^  m^s“^)  and  also  around  800  km  when 
Aft  is  reduced  to  8.0  x  10^°m‘^s“^.  For  equilibrium  C  the  jet  penetrates  1250  km  at 
both  its  bifurcation  point  (Aft  =  2.4  x  10^^  m^s“^)  and  at  A;,  =  8.0  x  10^”  m^s“^.  For 
equilibrium  D  the  jet  penetrates  up  to  1800  km  at  both  its  bifurcation  point  (Aft  = 
1.2  X  10^^  m^s~^)  and  at  A^  =  8.0  x  10^°  m^s“^  In  comparison,  the  anti-symmetric 
equilibrium  has  a  jet  and  recirculation  cells  whose  penetration  scale  increases  rapidly  as 
the  viscosity  coefficient  is  decreased.  Figure  3.25  shows  the  nearly  anti-symmetric  branch 
with  a  jet  which  penetrates  more  than  2000  km  into  the  interior  before  fanning  out  while 
the  jet  penetrated  only  400  km  at  the  first  pitchfork  bifurcation  PFa-  Consistent  with 
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this  behavior,  the  pairs  of  equilibria  which  bifurcated  at  the  highest  values  of  viscosity 
generally  have  a  jet  with  strong  meandering  but  with  only  weak  penetration  into  the 
interior.  On  the  other  hand,  the  solution  branch  E'  and  the  pairs  of  equilibria  which  have 
bifurcated  at  the  lower  values  of  viscosity  have  a  jet  with  little  meandering  but  which 
penetrates  strongly  into  the  interior  (see  Section  3.3.15).  All  the  solutions  are  identical  in 
the  region  away  from  the  recirculation  region  where  they  are  either  in  Sverdrup  balance, 
or  in  an  inertial  balance  with  the  streamlines  entering  the  western  boundary  layer. 

The  meandering  of  streamlines  across  the  line  of  zero  wind-stress  curl  (Figures  3.18 
through  3.25),  has  important  consequences  for  both  the  energy  and  vorticity  budgets  of 
the  steady  state  solutions,  which  we  now  address  in  more  detail. 

3.3.13  Global  Vorticity  Balance 

In  this  section  we  look  at  global  vorticity  budgets  for  each  gyre.  If  the  streamline  sep¬ 
arating  the  sub-polar  from  the  sub-tropical  gyre  is  not  coincident  with  the  line  of  zero 
wind-stress  curl,  the  circulation  can  advect  negative  vorticity  into  a  region  of  positive 
wind-stress  curl  and  vice-versa.  We  can  think  of  this  advection  of  vorticity  as  an  inter¬ 
gyre  flux  of  vorticity  provided  we  deflne  the  gyres  to  be  the  regions  occupied  by  the 
sub-tropical  and  sub-polar  gyres  of  the  linear  Munk-like  solution.  From  this  point  of 
view,  the  region  occupied  by  the  gyres  is  flxed,  and  consequently,  the  vorticity  input  by 
the  wind-stress  curl  is  also  flxed.  Tables  3.1  and  3.2  give  the  gyre  integrated  vorticity 
budget  for  the  sub-tropical  and  sub-polar  gyres  respectively.  The  advection  terms  in  the 
vorticity  equation  cannot  generate  any  vorticity;  they  act  to  only  redistribute  it.  Thus 
any  net  basin  integrated  input  of  vorticity  by  the  wind  must  be  removed  by  the  explicit 
friction  terms.  For  the  conflguration  used  in  this  chapter,  with  a  =  0.05,  the  net  input  of 
vorticity  by  the  wind  in  the  sub-tropical  gyre  is  -2.23  x  lO'^  s-^,  and  for  the  sub-polar 
gyre  it  is  2.0177  x  10“^  s“^.  The  sub-tropical  gyre  receives  5%  more  vorticity  from  the 
wind  than  the  sub-polar  gyre. 

Tables  3.1  and  3.2,  show  that  the  inter-gyre  flux  of  vorticity  is  crucial  for  equilibria 
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A,  B  and  B'  which  are  the  first  to  bifurcate.  Since  these  equilibria  are  furthest  in 
parameter  space  from  their  bifurcation  points,  they  are  the  least  anti-symmetric.  Also, 
the  solutions  B',  C,  and  Z)',  which  have  a  jet  that  first  turns  south  after  separating  from 
the  western  wall  have  weaker  inter-gyre  fluxes  of  vorticity  than  their  nearly  mirror  image 
counterparts  B,  C  and  D  which  have  a  jet  which  first  turns  north.  This  asymmetry  is 
due  to  the  weaker /stronger  vorticity  input  in  the  sub-polar/sub-tropical  gyre. 

The  dominant  explicit  dissipation  term  in  the  vorticity  equation  is  the  biharmonic 
viscosity.  It  generally  becomes  more  important  for  the  more  anti-symmetric  solution; 
although  the  relative  diflPerences  are  small  compared  with  the  changes  in  the  advection  and 
bottom  friction  terms.  The  sink  of  vorticity  through  lateral  diffusion  is  generally  stronger 
for  the  unprimed  solutions.  This,  along  with  the  weaker  inter-gyre  flux  of  vorticity 
for  the  unprimed  solutions  increases  the  importance  of  bottom  drag  for  removing  the 
excess  vorticity.  To  compensate  for  the  weaker  inter-gyre  flux  of  vorticity,  the  more  anti¬ 
symmetric  solutions  dissipate  much  more  vorticity  through  bottom  friction  than  do  the 
more  non-symmetric  solutions.  For  example,  bottom  friction  is  44%  more  important  for 
solution  E'  than  it  is  for  solution  C .  It  can  also  be  noticed  that  for  the  primed  solutions, 
bottom  friction  is  always  stronger  than  for  the  unprimed  counter  parts.  This  is  consistent 
with  the  weaker  inter-gyre  vorticity  flux,  and  weaker  lateral  diffusion.  Interestingly,  the 
time  dependent  trajectory  is  most  like  the  primed  solution,  with  a  quasi-permanent 
southward  meander  shortly  after  the  jet  separates  from  the  coast. 

3.3.14  Global  Energy  Balance 

In  addition  to  having  different  flow  fields,  each  equilibrium  state  has  a  different  energy 
level.  The  difference  in  the  energy  level  maintained  by  each  state  is  due  to  the  fact  that 
both  the  energy  dissipation  and  the  energy  input  by  the  wind  stress  are  a  function  of  the 
flow  field.  The  energy  input  by  the  wind-stress  is  given  by  the  correlation  between  the 
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Solution 

Vorticity  input 
(10-3s-2) 

Bottom  dissipation 
(l0-3s-2) 

Lateral  diffusion 
(10-"s-2) 

Advection 

(l0-3s-2) 

Linear 

-2.2300 

-0.5248 

-1.7052 

-0.0000 

A 

-2.2300 

-0.3411 

-0.8845 

-1.0044 

B 

-2.2300 

-0.3513 

-0.8603 

-1.0184 

B' 

-2.2300 

-0.3904 

-0.9431 

-0.8965 

C 

-2.2300 

-0.3138 

-0.9001 

-1.0161 

C 

-2.2300 

-0.4492 

-0.9838 

-0.7970 

D 

-2.2300 

-0.3270 

-0.9087 

-0.9943 

D' 

-2.2300 

-0.6193 

-0.9560 

-0.6547 

E' 

-2.2300 

-0.7051 

-0.9130 

-0.6120 

Table  3.1:  Integrated  vorticity  budget  for  sub-tropical  gyre. 


Solution 

Vorticity  input 
(10-V2) 

Bottom  dissipation 
(10-V2) 

Lateral  diffusion 
(10-3s-2) 

Advection 

(10-3s-2) 

Linear 

2.0177 

0.4757 

1.5420 

0.0000 

A 

2.0177 

0.2879 

0.7254 

1.0044 

B 

2.0177 

0.2934 

0.7058 

1.0184 

B' 

2.0177 

0.3405 

0.7807 

0.8965 

C 

2.0177 

0.2471 

0.7545 

1.0161 

C 

2.0177 

0.3997 

0.8209 

0.7970 

D 

2.0177 

0.2602 

0.7632 

0.9943 

D' 

2.0177 

0.5653 

0.7976 

0.6547 

E' 

2.0177 

0.6454 

0.7603 

0.6120 

Table  3.2:  Integrated  vorticity  budget  for  sub-polar  gyre. 
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can  in  fact  be  attributed  to  differences  in  the  potential  energy  while  bottom  dissipation 
is  proportional  to  the  kinetic  energy.  For  comparison,  Table  3.3  also  gives  the  energy 
balance  for  the  linearized  model.  For  this  solution,  lateral  diffusion  accounts  for  more 
than  half  the  energy  dissipation.  Since  there  are  no  inertial  effects  for  the  linearized 
model,  all  streamlines  pass  through  the  frictional  boundary  layer.  The  absence  of  iner¬ 
tial  effects  also  prevents  recirculation  cells  from  forming,  thereby  eliminating  important 
regions  where  bottom  friction  dissipates  energy. 

3.3.15  Zonal  Jet  Penetration 

The  most  striking  difference  between  the  multiple  equilibria  is  the  extent  of  the  jet 
penetration  into  the  basin  interior  and  the  amount  of  meandering  by  the  streamlines. 
The  issue  of  the  zonal  penetration  scale  of  mid-latitude  oceanic  jets  is  an  interesting 
problem  in  its  own  right.  Holland  and  Schmitz  (1985),  Greatbatch  (1988),  Marshall 
and  Marshall  (1992),  among  others  have  addressed  the  problem  of  what  sets  the  zonal 
penetration  scale  of  mid-latitude  jets  in  ocean  models.  Marshall  and  Marshall  as  well  as 
Greatbatch  argued  that  for  free  jets  with  q  =  F('0),  the  sign  of  dq/dip  controls  whether 
the  jet  strikes  seaward  or  turns  back  on  itself.  They  drew  an  analogy  between  two 
analytic  free  solutions  of  the  barotropic  vorticity  equation,  the  modon  and  the  Fofonoff 
mode,  and  the  penetration  characteristic  of  a  free  mid-latitude  jet.  Both  the  modon 
and  the  Fofonoff  mode  have  linear  ?  -  ^  relationships.  The  modon  with  dq/d^)  <  0 
has  a  circulation  which  is  confined  by  a  bounding  streamline  beyond  which  there  is  no 
flow.  The  Fofonoff  mode  on  the  other  hand,  has  dq/dxj)  >  0  and  a  circulation  that  fills 
the  entire  basin.  Based  on  these  analogies,  Marshall  and  Marshall  (1992),  went  on  to 
perform  a  series  of  numerical  experiments  in  which  they  controlled  the  9  -  -0  relationship 
for  the  jet  by  specifying  inflow  conditions,  and  showed  that  indeed,  when  the  prescribed 
inflow  conditions  had  dqjd’ip  <  0,  the  jet  turned  back  on  itself  and  when  dq/dx!)  >  0,  the 
jet  penetrated  across  the  basin  reaching  the  eastern  wall.  Based  on  these  experiments, 
they  proposed  the  theory  that  the  different  penetration  scale  for  the  mid-latitude  jets  in 
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different  models  can  be  explained  by  the  different  sign  of  dq/dil^. 

From  the  steady  state  solutions  presented  in  this  chapter,  it  is  however  possible  to 
state  that  the  theory  of  Marshall  and  Marshall  (1992)  is  incorrect  In  fact,  the  solutions 
found  here  have  jets  with  very  different  penetration  scale  while  having  nonetheless  similar 
q  —  tp  relationships.  For  free  solutions,  Vg  and  are  parallel  so  that  the  following 
relationship  holds 

|=sign(V,.V^)M  (3.18) 

The  sign  of  (Vg  •  V^)  can  be  computed  for  the  numerical  solutions  we  have  found,  and 
to  the  extent  that  Vq  and  Vtp  are  aligned,  it  gives  the  sign  of  dq/d'tp.  This  was  done  for 
each  of  the  solutions.  Figures  3.26(a),  3.27(a),  3.28(a)  and  3.29(a)  show  the  results  for 
equilibria  A,  B',C',  and  E'.  Shaded  areas  show  regions  in  which  dq/dip  >  0.  Also  shown 
are  the  ^-contours  for  the  solution.  Figures  3.26(b),  3.27(b),  3.28(b)  and  3.29(b)  show 
regions  in  which  the  g-contours  are  nearly  parallel  to  the  h  contours.  Shaded  areas  re 
regions  in  which  the  angle  between  Vq  and  Vh  is  less  than  10®.  Also  shown  are  the  h 
contours.  We  can  see  that  the  recirculation  cells  always  have  dq/dip  positive  regardless 
of  the  penetration  extent  of  the  jet.  This  is  in  contradiction  to  the  theory  of  Marshall 
and  Marshall  (1992),  which  predicts  a  negative  dq/dip  within  the  recirculation  cells  for 
jets  which  turn  back  on  themselves. 

The  sign  of  Vq  •  V‘ip  plotted  in  Figures  3.26  through  3.29  can  be  used  to  understand 
why  some  regions  have  stationary  waves,  and  others  do  not.  Recall  that  Vip  gives  the 
direction  in  which  the  current  flows  and  that  Vq  gives  the  direction  in  which  Rossby  waves 
propagate  their  phases.  Regions  which  have  Vq  •  Vip  negative  can  support  stationary 
waves  because  the  flow  opposes  the  waves.  Consistent  with  this.  Figures  3.26  through  3.29 
show  that  regions  with  Vq  ■  Vip  <  0  (unshaded)  have  wavy  g-contours.  The  waviness  of 
the  solutions  to  the  east  of  the  recirculation  cells  is  qualitatively  similar  to  the  dynamical 
regime  defined  by  Moore’s  solution. 
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Figure  3.26:  (a)  Contour  plot  of  potential  vorticity  (q)  for  equilibrium  A.  Shaded  areas 
are  regions  in  which  Vq-Vip  >  0,  and  unshaded  areas  are  regions  where  Vq  •  <  0. 

(b)  Contour  plot  of  the  corresponding  interface  height  anomaly  (h).  Shaded  areas  are 
regions  in  which  the  angle  between  V?  and  Vh  is  less  than  10°. 


112 


Figure  3.27:  (a)  Contour  plot  of  potential  vorticity  for  equilibrium  B  .  Shaded  areas 
are  regions  in  which  Vq  ■  >  0,  and  unshaded  areas  are  regions  where  Vq'  •  "Vip  <  0. 

(b)  Contour  plot  of  the  corresponding  interface  height  anomaly  (h).  Shaded  areas  are 
regions  in  which  the  angle  between  V?  and  Vh  is  less  than  10°. 
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Figure  3.28:  (a)  Contour  plot  of  potential  vorticity  for  equilibrium  C.  Shaded  areas  are 
regions  in  which  Vq  ■  Vtp  >  0,  and  unshaded  areas  are  regions  where  Vg  •  Vijj  <  0.  (b) 
Contour  plot  of  the  corresponding  interface  height  anomaly  (h).  Shaded  areas  are  regions 
in  which  the  angle  between  Vg  and  V/i  is  less  than  10”. 
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3.4  Fixed  Points  and  Time-dependent  simulations 

In  Section  3.2  we  described  the  phenomenology  of  the  time  dependent  simulation.  In 
this  section  we  will  revisit  the  time-dependent  behavior  of  the  model  and  explore  the 
possibility  that  there  is  some  connection  between  the  trajectory  of  the  model  in  phase 
space  and  the  model’s  fixed  points. 

There  is  a  remarkable  similarity  between  the  contour  plots  of  the  interface  anomaly 
for  the  flow  averaged  within  the  high  energy  regime  and  the  fixed  points  E'  (compare 
Figure  3.25  to  Figure  3.5),  for  the  medium  energy  regime  and  the  fixed  point  C  (compare 
Figure  3.22  and  Figure  3.5)  and  for  the  low  energy  regime  and  the  fixed  point  B'  (compare 
Figure  3.20  and  Figure  3.5).  The  zonal  penetration  scale  for  the  jet  is  roughly  equal  for 
the  fixed  points  and  for  the  time  averaged  flows.  Also,  the  meandering  structure  of  the 
jet  is  remarkably  similar. 

It  is  useful  to  first  compare  the  energy  level  of  the  time-dependent  simulation  and 
the  energy  levels  of  the  fixed  points.  Figure  3.30  shows  a  histogram  of  the  total  energy 
for  the  1500  year  time  series  as  well  as  the  energy  level  for  each  steady  state  solution. 
Very  little  of  the  distribution  density  spreads  to  energy  levels  higher  than  the  level  of 
equilibrium  E'.  As  we  will  show  below  most  of  the  density  of  high  energy  realizations  can 
be  attributed  to  trajectories  which  tend  towards  equilibrium  E'  from  low  energy  levels. 

Except  for  the  general  agreement  between  the  order  of  magnitude  of  the  energy  levels 
of  the  fixed  points  and  the  time  dependent  trajectory,  there  is  no  clear  agreement  between 
the  energy  levels  of  the  steady  state  solutions  and  the  peaks  in  the  histogram  near 
3.55  X  10"^^  J  (low  energy)  and  3.95  x  10“^'’'  J  (medium  energy).  Equilibrium  C  is 
close  to  the  peak  in  the  distribution  near  3.95  x  10^^  J,  but,  as  we  shall  see  using  the 
distance  diagnostic  in  Section  3.4.4,  when  the  model  trajectory  is  closer  to  equilibrium 
C  its  energy  level  is  more  often  in  the  range  of  the  low  energy  peak.  Similarly,  even 
though  it  appears  that  equilibrium  D  and  D'  line  up  with  a  possible  high  energy  peak 
near  4.3  x  10  J,  when  the  model  is  in  this  energy  range  it  is  usually  much  closer  to 
equilibrium  E'.  Finally,  when  the  model  is  closest  to  equilibrium  D',  it  is  usually  in  the 
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Figure  3.31:  Time  averaged  interface  height  anomaly  field  (C.I.  20  m). 

medium  or  low  energy  level,  despite  the  fact  that  the  energy  level  of  this  fixed  point  is 
in  the  range  of  the  high  energy  regime. 

3.4.1  Fixed  points  and  Modes  of  Variability 

If  the  idea  that  the  model’s  fixed  points  act  to  steer  the  time-dependent  trajectory  in 
phase  space  is  correct,  we  would  expect  to  see  modes  of  variability  associated  with  struc¬ 
tures  in  phase  space  that  point  away  from  the  time-mean  state  and  towards  the  fixed 
points.  We  will  restrict  the  study  to  the  fixed  points  E',  D',  C  and  B'  since  they  are  most 
similar  to  the  time  averaged  flows  computed  by  McCalpin  and  Haidvogel.  In  Figure  3.31 
the  time-mean  interface  height  anomaly,  h,  is  contoured.  It  is  obtained  by  averaging  the 
field  saved  at  5  day  intervals  over  a  period  of  1200  years  excluding  the  spin-up  period. 
The  amount  of  variability  away  from  this  mean  state  and  towards  the  fixed  points  E', 
D  ,  C  ,  and  B  ,  can  be  evaluated.  Using  the  Gram-Schmidt  orthonormalization  process, 
four  orthonormal  vectors  which  span  the  directions  in  phase  space  that  point  away  from 
the  mean  state  and  towards  the  four  primed  fixed  points  can  be  created.  The  first  mode. 
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model  in  Figure  3.32,  is  simply  the  normalized  difference  between  the  mean  state  and 
fixed  point  E': 

vV  =  h-E' 
model  =  vl' / {vV ,  vl') . 

The  second  mode,  mode2  in  Figure  3.32,  is  the  normalized  difference  between  the  mean 
state  and  fixed  point  D',  less  the  projection  onto  model: 

v2'  =  h-D' 
v2"  =  v2' —  (Ji,  model) 
mode2  =  v2"  I  {v2"  ,v2"). 

The  third  mode,  mode?)  in  Figure  3.32,  is  the  normalized  difference  between  the  mean 
state  and  fixed  point  C",  less  the  projections  onto  model  and  mode2: 

V?'  =  h-C 

V?"  =  V?'  —  (h,  model)  —  (h,  mode2) 
mode?  =  V?"  /  {v?"  ,v?"). 

Finally  the  fourth  mode,  modeA  in  Figure  3.32,  is  the  normalized  difference  between  the 
mean  state  and  fixed  point  B',  less  the  projection  onto  the  three  previous  modes: 

vA'  =  h-B' 

vA"  =  vA'  —  (h,  model)  —  (h,  mode2)  —  (h,  mode?) 
modeA  =  t;4"). 

In  Figure  3.32  we  show  the  contour  plots  of  model,  mode2,  mode?  and  modeA.  The 
amount  of  variability  associated  with  each  mode  is  obtained  by  projecting  the  deviations 
away  from  the  mean  onto  the  orthonormal  set.  Approximately  30  %  of  the  total  variability 
is  captured  by  the  four  modes.  The  amount  is  significant  considering  that  the  system  as 
24881  degrees  of  freedom.  The  break-up  of  the  variability  associated  with  each  mode  is 
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Figure  3.32:  Four  orthonormal  modes  (model,  mode2,  modeS,  and  modeA)  spanning  the 
hyper  plane  defined  by  the  difference  from  the  time-averaged  flow  field  and  the  steady 
state  solutions  E',  D' ,  C",  and  B' 


as  follows:  model,  which  points  in  the  direction  of  E'  captured  11%  of  the  variability, 
mode2  captured  7.9%  of  the  variability,  mode3  captured  7.2%  of  the  variability,  and 
modeA  captured  4.4%  of  the  variability.  Furthermore,  most  of  the  variance  in  the  interface 
height  anomaly  capture  by  the  four  modes,  is  at  low  frequencies.  Figure  3.33  shows  a 
plot  of  the  frequency  times  power  density  spectrum  for  the  basin  integrated  variance  of 
the  interface  height  anomaly  for  the  full  field,  and  for  the  field  in  which  the  projection 
onto  the  four  modes  has  been  removed.  The  plot  shows  a  significant  part  of  the  variance 
associated  with  periods  longer  than  1  year  project  onto  the  four  modes.  The  fixed  points 
of  the  model  do  appear  to  steer  the  model  trajectory  in  phase  space,  and  thereby  generate 
modes  of  low-frequency  variability. 
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Interface  Height  Anomaly  Variance  Spectrum 


Figure  3.33;  Frequency  times  power  density  spectrum  for  the  variance  of  the  interface 
height  anomaly  of  the  total  field  (upper  curve),  and  for  the  field  in  which  the  part  of 
the  variance  which  projects  onto  model,  mode2,  modeS,  and  modei  has  been  removed 
(lower  curve).  The  dashed  lines  are  95%  confidence  intervals. 
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3.4.2  Distance  Diagnostic 

To  quantify  in  a  more  objective  manner  the  similarity  between  the  fixed  points  and  the 
time  averaged  flows  within  each  of  the  high,  medium  and  low  energy  regimes,  the  distance 
in  phase  space  between  a  fixed  point  and  the  instantaneous  model  state  can  be  evaluated. 
The  distance  dx{t)  between  a  fixed  point  X  and  the  model  state  at  time  t  is  give  by 

4'W  =  \po9'  J f  -  X{x,y)fdxdy,  (3.19) 

in  which  potential  energy  was  used  as  the  norm. 

Using  the  1500  year  time-series  of  the  interface  height  anomaly,  sampled  at  5  day 
intervals,  8  time-series  for  the  distance  between  the  instantaneous  model  state  and  the  8 
fixed  points  were  computed.  The  distance  diagnostic  reveals  a  more  complicated  picture 
than  the  simple  interpretation  that  there  is,  for  example,  a  one  to  one  correspondence 
between  the  high,  medium,  and  low  energy  regimes  and  the  fixed  points  E',  C  and  B'. 
Before  discussing  the  results,  we  first  illustrate  the  use  of  the  distance  diagnostic  with 
the  familiar  Lorenz  model. 

3.4.3  Example:  Distance  Diagnostic  for  the  Lorenz  Model 

In  this  subsection  we  use  the  Lorenz  model  to  illustrate  the  idea  of  distance  in  phase 

the  following  set  of  three  ordinary  differential 

(r{Y-X),  (3.20) 

rX-Y-  XZ,  (3.21) 

XY  -  bZ.  (3.22) 

With  the  parameters  chosen  as  follows:  cr  =  10,  r  =  28,  and  b  =  8/3,  the  model 
trajectories  approach  the  well  known  Lorenz  attractor  (Figure  3.34).  In  Figure  3.34  the 
circles  indicate  the  fixed  points  of  the  model  which  are  located  at  {X  =  0,Y  =  0,Z  =  0) 


space.  The  Lorenz  model  is  given  by 
equations, 

dX 

dt 

d^ 

dt 

dZ 

dt 
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Figure  3.34;  Typical  trajectory  of  the  Lorenz  model  illustrating  the  Lorenz  attractor. The 
three  circles  indicate  the  three  fixed  points  ((7+, (7“  and  O). 

labeled  O,  at  {X  =  6\/2,  Y  =  6\/2,  Z  =  27)  labeled  (7+,  and  at  {X  =  -6\/2,  Y  =  -6V2, 
Z  =  27)  labeled  C~.  The  model  attractor  consists  roughly  of  two  spiraling  lobes  centered 
on  the  fixed  points  C"*"  and  C~ .  The  model  trajectory  exhibits  irregular  transitions  from 
one  lobe  to  the  other. 

Figure  3.35  shows  a  typical  time  series  of  dc+{t),  the  distance  from  fixed  points  C~^ 
to  the  model  trajectory  and  dc-{t),  the  distance  from  fixed  points  C~  to  the  model 
trajectory.  In  general,  when  dc+  is  small,  dc-  is  relatively  larger  and  executes  large 
amplitude  oscillations.  Similarly,  when  when  dc-  is  small,  it  is  dc+  which  is  larger,  and 
which  executes  large  amplitude  oscillations. 

The  distance  diagnostic  can  be  used  to  illustrate  the  property  that  the  closer  the 
model  trajectory  gets  to  one  of  the  fixed  points  or  C~,  the  longer  the  trajectory 
will  spiral  in  the  lobe  centered  on  the  fixed  point.  Figure  3.36  shows  a  scatter  plot  of 
the  minimum  distance  from  the  instantaneous  model  state  to  the  fixed  point  achieved 
for  each  event  where  the  model  spirals  on  the  lobe  centered  on  versus  the  duration 
of  the  event.  We  define  an  event  to  be  an  uninterupted  period  of  time  in  which  the 
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Figure  3.35:  (a)  Typical  time  series  of  the  distance  from  the  model  trajectory  to  the  fixed 
point  (7+.  (b)  Typical  time  series  of  the  distance  from  the  model  trajectory  to  the  fixed 
point  C~ . 

instantaneous  model  trajectory  is  closer  to  one  fixed  point  rather  than  the  other.  The 
duration  of  each  event  increases  as  the  minimum  distance  to  the  fixed  point  decreases. 
As  will  be  shown  in  the  next  section,  this  is  a  feature  shared  by  one  of  the  fixed  points 
of  the  ocean  model. 

The  fixed  points  of  the  Lorenz  model  can  be  thought  of  as  analogs  of  the  “basic 
state  of  the  ocean  model  within  the  high  medium  and  low  energy  regimes,  and  the 
spirals  around  the  fixed  points  of  the  Lorenz  model  are  analogs  of  the  eddy  variability 
associated  with  each  regime.  In  the  next  section  we  will  test  whether  this  analogy  can 
be  made  more  concrete  by  comparing  the  fixed  points  of  the  ocean  model  to  the  model’s 
attractor. 

3.4.4  Distance  Diagnostic  Applied  to  the  Ocean  Model 

If  we  define  the  model  trajectory  to  be  in  regime  Rx  whenever  it  is  closer  to  fixed  point 
X  than  to  any  other  fixed  point,  it  is  possible  to  compute  the  fraction  of  time  in  which 
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Figure  3.36:  Scatter  plot  of  the  minimum  distance  to  the  fixed  point  as  a  function  of 
the  duration  time  for  which  the  model  trajectory  is  in  the  lobe  centered  on  . 

the  model  state  will  be  simultaneously  in  regime  Rx  and  the  low  {TE  <  3.75  x  lO^^J), 
medium  (3.75  x  10^^ J  <TE  <  4.25  x  lO^’^J)  and  high  {TE  >  4.25  x  lO^’^J)  energy  states 
defined  by  McCalpin  and  Haidvogel  (1996). 

Table  3.4  gives  the  fraction  of  the  time  that  the  model  is  in  a  given  regime  (defined 
by  its  proximity  to  a  steady  state)  and  is  also  in  either  the  high,  medium  or  low  energy 
states.  In  general,  the  fixed  point  E'  is  associated  with  high  and  medium  energy  levels, 
the  fixed  point  D'  is  associated  with  low  and  medium  energy  levels,  the  fixed  point  C 
is  associated  with  medium  and  low  energy  levels  and  the  fixed  point  B'  is  associated 
only  with  low  energy  levels.  The  table  shows  that  when  the  model  is  in  the  high  energy 
range  it  is  usually  in  regime  Re'-  There  is  however  a  nonzero  probability  that  it  will 
be  in  either  regime  Rd  or  Re'-  When  the  model  is  in  the  Medium  energy  regime  it  is 
usually' in  regimes  Rd'  and  Re',  and  when  the  model  is  in  the  low  energy  range,  there  is 
a  non-negligible  probability  that  it  will  at  some  time  be  closer  to  any  given  fixed  point 
than  to  any  other,  although  it  is  usually  closest  to  fixed  points  C  and  B' . 

If  the  cross  tabulations  are  restricted  to  persistent  regimes,  a  simpler  picture  emerges 
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Regime 

High 

Medium 

Low 

Ra 

0.00 

0.01 

0.07 

Rb 

0.00 

0.01 

0.03 

Rb' 

0.00 

0.00 

0.22 

Rc 

0.00 

0.01 

0.01 

Re 

0.00 

0.16 

0.35 

Rd 

0.09 

0.04 

0.03 

Rd' 

0.03 

0.47 

0.18 

Re' 

0.88 

0.31 

0.10 

Table  3.4:  Fraction  of  time  spent  closest  to  each  of  the  fixed  points  given  that  it  is  either 
in  the  High,  the  Medium  or  the  Low  Energy  regime. 


Regime 

3-5  years 

5-10  years 

>  10 

Ra 

0 

0 

0 

Rb 

0 

0 

0 

Rb' 

0 

0 

0 

Rc 

0 

0 

0 

Rc 

1 

0 

0 

Rd 

0 

0 

0 

Rd' 

21 

1 

0 

Re' 

3 

11 

9 

Table  3.5:  Number  of  events  in  which  the  model  stayed  in  a  particular  regime  for  periods 
of  time  between  3  and  5  years,  between  5  and  10  years  and  more  than  10  years. 


for  the  role  of  fixed  point  E' .  In  Table  3.5  we  give  the  number  of  realizations  of  a  given 
regime  as  defined  by  the  proximity  to  a  fixed  point  for  a  duration  of  3-5  years,  5-10  years, 
and  longer  than  10  years.  From  this  table  we  can  see  that  only  regime  Re>  persists  for 
periods  of  time  greater  than  3  years.  For  comparison.  Table  3.6  shows  the  number  of 
occurrences  of  the  high,  medium  and  low  energy  regimes  that  persist  for  lengths  of  time 
between  3-5  years,  5-10  years,  and  longer  than  10  years.  There  is  an  exact  correspondence 
between  events  where  the  model  state  is  in  the  high  energy  regime,  and  that  the  model 
is  closest  to  equilibria  E' .  However,  when  the  model  trajectory  is  in  the  medium  and  low 
energy  regimes,  for  persistent  periods  of  time,  it  does  not  stay  persistently  in  the  regimes 
associated  with  any  particular  fixed  point. 

The  correspondence  between  the  various  regimes  as  defined  by  the  model’s  energy 
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Regime 

<3  —  5  years 

5—10  years 

>  10 

Low 

23 

14 

9 

Medium 

23 

44 

33 

High 

3 

11 

9 

Table  3.6:  Number  of  events  in  which  the  model  stayed  in  a  particular  regime  for  periods 
of  time  between  3  and  10  years,  between  5  and  10  years,  and  more  than  10  years. 


Regime 

<3  —  10  years 

10  —  15  years 

>15  years 

Ra 

0.1400 

0.0912 

0.0211 

Rb 

0.0415 

0.0738 

0.0426 

Rb' 

0.2385 

0.3529 

0.2965 

Rc 

0.0369 

0.0000 

0.0000 

Re 

0.3908 

0.3442 

0.5017 

Rd 

0.0600 

0.0000 

0.0000 

Rd' 

0.0908 

0.1086 

0.1242 

Re' 

0.0015 

0.0293 

0.0110 

Table  3.7:  Fraction  of  the  time  spent  closer  to  a  given  fixed  point  than  to  any  other  fixed 
point,  given  the  the  model  is  in  the  Low  energy  regime  for  an  extended  period  of  time 

level  and  the  proximity  of  the  model  trajectory  to  the  various  fixed  points  can  be  further 
characterized.  For  persistent  regimes,  Table  3.7,  Table  3.8,  Table  3.9,  give  the  relative 
time  spent  closer  to  any  given  fixed  point  while  also  being  in  either  in  the  low,  medium 
or  high  energy  regimes.  From  Table  3.7  we  see  that  whenever  the  model  is  in  the  low 
energy  regime  for  persistent  periods  of  time,  the  model  trajectory  partitions  its  time  in 
the  regimes  Rci  Rb'  Rd'  and  Ra-  From  Table  3.8,  we  see  that  when  the  model  state  is 
in  the  medium  energy  regime,  it  partitions  its  time  in  the  Re,  Rd',  and  Re'  regimes.  For 
those  events  when  the  model  stays  in  the  medium  energy  level  for  more  than  10  years, 
it  spends  a  considerable  amount  of  time  in  regime  Re'-  Finally  from  Table  3.9,  we  see  a 
clear  connection  between  persistent  high  energy  events,  and  the  proximity  of  the  model 
trajectory  to  the  fixed  point  E'. 

In  Figure  3.37  we  plot  the  square  of  the  minimum  distance  to  the  fixed  point  E'  for 
the  23  persistent  high  energy  events  listed  in  Table  3.6  as  a  function  of  the  duration  time 
of  the  corresponding  events.  From  the  plot  we  see  that  the  closer  the  model  trajectory 
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Regime 

<3  —  10  years 

10  —  15  years 

>15  years 

Ra 

0.0223 

0.0000 

0.0000 

Rb 

0.0000 

0.0000 

0.0000 

Rb' 

0.0143 

0.0000 

0.0000 

Rc 

0.0048 

0.0000 

0.0000 

Re 

0.1927 

0.0249 

0.0112 

Rd 

0.0000 

0.0000 

0.0000 

Rd' 

0.6226 

0.5748 

0.3844 

Re> 

0.1433 

0.4002 

0.6045 

Table  3.8:  Fraction  of  the  time  spent  closer  to  a  given  fixed  point  than  to  any  other  fixed 
point,  given  the  the  model  is  in  the  Medium  energy  regime  for  an  extended  period  of 
time 


Regime 

<3  —  10  years 

10  —  15  years 

>15  years 

Ra 

0.0000 

0.0000 

0.0000 

Rb 

0.0000 

0.0000 

0.0000 

Rb' 

0.0000 

0.0000 

0.0000 

Rc 

0.0000 

0.0000 

0.0000 

rc 

0.0000 

0.0000 

0.0000 

Rd 

0.0000 

0.0000 

0.0000 

Rd' 

0.0000 

0.0000 

0.0000 

Re' 

1.0000 

1.0000 

1.0000 

Table  3.9:  Fraction  of  the  time  spent  closer  to  a  given  fixed  point  than  to  any  other  fixed 
point,  given  the  the  model  is  in  the  high  energy  regime  for  an  extended  period  of  time 
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Figure  3.37:  Scatter  plot  of  the  square  of  the  minimum  distance  to  the  fixed  point  E'  as 
a  function  of  the  duration  time  for  of  each  of  the  corresponding  high  energy  events  listed 
in  Table  3.6. 

gets  to  the  fixed  point  the  longer  the  high  energy  event  persists.  This  suggest  that  a 
more  detailed  study  of  the  structure  of  phase  space  near  the  steady  state  E'  might  help 
further  understand  the  dynamics  of  the  high  energy  events. 

3.5  Structure  of  Phase  Space  near  E' 

We  have  seen  in  the  previous  section  that  part  of  the  model  attractor  lies  close  the  fixed 
point  E'  in  phase  space.  Because  of  the  quasi-stability  of  the  high  energy  regime,  and  the 
clear  correspondence  between  persistent  high  energy  regimes  and  the  steady  state  solution 
E',  it  is  worthwhile  to  further  explore  the  behavior  of  the  model  in  the  neighborhood  of 
this  fixed  point. 
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3.5.1  Linear  Stability  Analysis 

A  linear  stability  analysis  of  the  fixed  point  E'  was  performed  following  the  method 
outlined  in  Appendix  B.  This  involves  looking  for  modal  solutions,  of  the  form 

y{x,  y,  t)  =  M{x,  y)e‘^\  (3.23) 

where  M (x,  y)  and  a  can  be  complex.  The  mode  can  be  rewritten  in  the  following  form 

V  =  Re{A{x,  =  A{x,  y)e^^cos{(j){x,  y)  +  ut),  (3.24) 

where  uj,  7,  <f>[x,y)  and  A{x,y)  are  all  purely  real.  The  quantity  7  is  the  growth-rate  of 
the  mode.  If  it  is  positive,  the  solution  is  unstable  and  the  mode  will  grow.  The  field 
A{x,y)  is  positive  everywhere,  and  gives  the  spatial  envelope  or  amplitude  of  the  mode. 
The  field  4>{x,y)  varies  from  -n/2  to  +7r/2  and  forms  part  of  the  phase,  (l){x,y)  -  ujt. 
Locally  we  can  approximate  (j){x,  y)  by  its  Taylor  expansion 

y)-(t)Q  +  =  (l)o  +  kx  +  ly,  (3.25) 

so  that  gives  the  direction  of  phase  propagation,  and  c{x,y)  =  uj/\V(f){x,y)\  gives 
the  local  phase  speed.  The  stronger  the  gradient  in  (j),  the  slower  the  phase  speed  and 
conversely,  the  smaller  the  gradient  in  (j)  the  faster  the  phase  speed. 

The  stability  analysis  reveals  that  the  fixed  point  E'  is  unstable  to  a  single  oscillating 
mode  with  a  period  of  approximately  T  =  21^ fu)  =  670  days,  and  a  growth-rate  with 
an  e-folding  time  of  I/7  =  897  days.  Figure  3.38(a)  shows  the  amplitude,  A{x,y),  of 
the  unstable  eigen-mode,  and  Figure  3.38(b)  shows  the  spatial  structure,  ^(rr,?/),  of  the 
unstable  eigen-mode’s  phase.  The  colormap  in  Figure  3.38(b)  varies  from  black  to  light 
gray  to  dark  gray  as  the  phase  varies  from  -7r/2  to  7r/2,  i.e  through  one  half  cycles. 

The  amplitude  of  the  unstable  mode,  (Figure  3.38(a)),  has  a  spatial  structure  con¬ 
sisting  of  4  cells  located  along  the  jet  axis  and  weaker  cells  situated  over  the  recirculating 
currents  to  the  north  and  south.  The  strongest  cell  is  situated  near  1700  km  eastward 
of  the  western  wall.  The  next  strongest  cell  is  situated  at  1150  km  and  its  intensity  is 
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Amplitude:  A(x,y) 


Phase:  <{>(x,y) 


X(km) 

Figure  3.38:  (a)  Amplitude  of  the  unstable  eigen-mode,  (b)  Phase  of  the  unstable  eigen¬ 
mode.  The  phase  propagates  from  black,  to  light  gray  to  to  dark  gray. 


60%  that  of  the  strongest.  The  next  cell  is  centered  at  600  km  and  its  intensity  is  60% 
that  of  the  strongest.  The  fourth  cell,  situated  close  to  the  western  wall,  has  an  intensity 
which  is  50%  that  of  the  strongest  cell.  The  cells  situated  over  the  westward  flowing 
recirculation  currents  have  maximums  near  1900  km  east  of  the  western  wall  and  decay 
slowly  to  the  west  and  rapidly  to  the  east.  The  strength  of  these  cells  is  approximately 
10%  of  the  strongest  cell  located  on  the  jet  axis. 

The  phase  structure,  (Figure  3.38(b))  of  the  unstable  mode  is  more  complicated.  Over 
most  of  the  basin,  where  the  basic  state  is  in  Sverdrup  balance,  and  where  the  amplitude 
of  the  mode  is  negligible,  the  phase  propagates  slowly  to  the  west,  wrapping  around 
the  jet  and  recirculation  gyres,  so  that  the  lines  of  constant  phase  become  progressively 
more  aligned  in  the  east-west  direction,  and  progress  towards  the  jet  from  the  north  and 
from  the  south.  In  the  region  of  the  westward  flowing  recirculating  currents,  the  phase 
progresses  westward  at  a  nearly  constant  speed,  (Figure  3.39(a), (c))  going  through  one 
full  cycle  as  it  progresses  from  the  eastern  tip  of  the  jet  to  the  western  boundary  current 
in  approximately  670  days.  The  phase  speed  of  approximately  0.045  ms~^  is  consistent 
with  the  phase  speed  of  long  Rossby  waves,  =  0.045  ms“\  and  with  the  westward 

drift  speed  of  rings  in  the  model.  McCalpin  and  Haidvogel  estimated  a  drift  speed  for  the 
rings  which  was  typically  within  10%  the  phase  speed  for  long  Rossby  waves.  Along  the 
jet  axis,  the  direction  of  phase  propagation  changes  directions  several  times  as  can  be  seen 
from  Figure  3.39(b),  and  there  are  several  extended  regions  where  the  phase  speed  nearly 
vanishes.  The  propagation  of  the  phases  with  varying  speed  under  the  modulation  of  the 
spatially  varying  amplitude  results  in  a  complicated  pattern  of  expanding  and  contracting 
cells  which  appear  to  be  stationary  at  times,  and  to  move  very  quickly  at  others. 

An  initial  value  problem  starting  from  fixed  point  E'  was  solved  numerically  for  20 
years.  Figure  3.40(a)  shows  the  evolution  of  the  potential  energy,  and  Figure  3.40(b) 
shows  the  evolution  of  the  kinetic  energy  for  the  first  20  years.  The  very  weak  growth-rate 
of  the  unstable  mode  explains  the  relative  stability  of  the  high  energy  regime.  Within 
a  span  of  20  years  to  total  energy  varies  by  less  than  1  %.  The  time- varying  structures 
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Figure  3.39:  (a)  Zonal  phase  speed  along  V  =  1600  km,  the  position  of  the  northern 
recirculation  current.  (b)Zonal  phase  speed  along  Y  =  1400  km,  the  position  of  the  jet. 
(c)  Zonal  phase  speed  along  V  —  1200  km,  the  position  of  the  southern  recirculation 
current.  The  speed  is  positive  eastward  and  negative  westward. 


Potential  Energy:  ic  =  E 


X  10’®  Kinetic  Energy:  ic  =  E 


Figure  3.40:  (a)  Time  series  for  the  potential  energy  for  a  simulation  starting  with  fixed 
point  E'  as  initial  condition,  (b)  Corresponding  time  series  for  the  kinetic  energy. 

observed  during  the  20  year  simulation  are  very  similar  to  those  of  the  unstable  eigen¬ 
mode,  and  to  the  structures  observed  in  full  simulation  when  the  the  model  trajectory  is 
near  the  fixed  point  E' . 

3.6  Discussion 

In  this  Chapter  we  have  demonstrated  that  much  of  the  low-frequency  variability,  30  %  of 
the  total  variance,  in  the  interface  height  anomaly  is  associated  with  structures  that  lie  in 
a  4-dimensional  hyper-plane  spanned  by  modes  formed  by  taking  the  difference  between 
the  mean  state  of  the  model  and  the  fixed  points  E',  D' ,  C  and  B'.  The  last  three  of 
these  equilibria  are  members  of  three  pairs  of  equilibria  which  are  the  result  of  three 
successive  symmetry  breaking  pitchfork  bifurcations.  The  asymmetry  in  the  wind-stress 
curl  has  selected  one  member  from  each  pair  of  nearly  mirror  image  equilibria,  namely, 
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the  one  with  a  jet  which  first  turns  south,  after  separating  from  the  western  wall  near 
the  line  of  zero  wind-stress  curl. 

Furthermore,  we  have  demonstrated  that  the  high  energy  regime  identified  by  Mc- 
Calpin  and  Haidvogel  (1996)  is  associated  with  the  existence  of  a  nearly  anti-symmetric 
fixed  point  or  steady  state  solution  which  we  labeled  E'.  A  modal  stability  analysis 
has  revealed  that  this  fixed  point  is  unstable  to  a  single  mode  which  is  oscillatory  in 
nature.  The  growth-rate  of  the  mode  is  weak,  with  an  e-folding  time  of  nearly  2.5  years, 
explaining  in  part  the  quasi-stability  of  the  high-energy  regime. 

McCalpin  and  Haidvogel  (1996),  have  noted  that  as  the  wind-stress  curl  profile  is 
made  progressively  more  asymmetric  by  increasing  a,  the  occurrence  of  persistent  the 
high  energy  events  decreases  gradually.  Consistent  with  this  observation,  the  existence 
of  the  fixed  point  E'  depends  on  the  parameter  a.  Only  for  —0.166  <  a  <  -t-0.166  does 
the  fixed  point  exist.  Figure  3.41  shows  a  bifurcation  diagram  of  the  total  energy  as  a 
function  of  a  for  the  fixed  points  E'  and  D'.  At  a  =  0.166  there  is  a  turning  point  where 
the  branch  E'  merges  with  the  branch  D'.  The  existence  of  this  turning  point  is  reflected 
in  the  time-dependent  simulations  of  McCalpin  and  Haidvogel.  They  found  that  by  the 
time  a  was  increased  to  0.1,  there  were  no  high  energy  events  which  persisted  for  periods 
greater  than  3  years. 

Finally,  the  important  role  played  by  the  basin  geometry  in  modifying  the  bifurcation 
structure  of  the  model  should  be  pointed  out.  The  larger  zonal  to  meridional  aspect 
ratio  of  the  basin  used  in  this  chapter  has  allowed  twice  as  many  pitchfork  bifurcations 
to  occur  occur  than  for  a  narrower  basin  (and  perhaps  other  pitchfork  bifurcations  would 
occur  if  the  biharmonic  viscosity  had  been  reduced  even  further).  In  addition,  the  wider 
basin  has  clarified  the  nature  of  the  symmetry  breaking  bifurcations,  by  making  evi¬ 
dent  the  quantized  nature  of  the  meanders  in  the  jet  separating  the  recirculation  gyres. 
Each  bifurcation,  adds  an  additional  half  meander  to  the  non-symmetric  solutions.  In 
the  calculation  by  Cessi  and  lerley  (1995)  the  narrow  basin  allows  only  two  pitchfork 
bifurcations  to  occur.  They  form  the  first  two  in  a  sequence  that  would  get  progressively 
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Figure  3.41:  Bifurcation  diagram  of  total  energy  in  Joules  versus  a.  There  is 
point  where  the  branches  for  fixed  points  E'  and  D'  merge  at  a  =  0.0166. 


a  turning 


longer  if  the  basin  was  allowed  to  get  progressively  wider.  Also,  the  calculations  presented 
here  suggest  that  the  interaction  of  the  jet  with  the  eastern  wall  might  be  responsible  for 
the  existence  of  the  saddle  node  bifurcation  leading  to  multiple  anti-symmetric  equilibria 
found  by  lerley  and  Sheremet  (lerley  and  Sheremet  1995)  for  the  single  gyre  case.  Again 
the  narrow  aspect  ratio  of  the  basin  obscures  this  possible  connection. 
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Chapter  4 


Conclusions 

4.1  Summary  of  the  Thesis 

The  major  hypothesis  presented  in  this  thesis  is  that  both  stable  and  unstable  steady  solu¬ 
tions  are  useful  in  describing  and  explaining  the  time-mean  and  low-frequency  variability 
of  ocean  models.  Steady  and  time-dependent  solutions  are  examined  for  two  ocean  mod¬ 
els  of  the  wind-driven  circulation.  The  first  model,  studied  in  Chapter  2,  is  a  barotropic 
model  driven  by  a  sinusoidal  wind-stress  curl  profile  in  a  rectangular  basin  with  a  zonal 
to  meridional  aspect  ratio  of  1:2.  The  explicit  dissipation  consists  of  bottom  drag  and 
biharmonic  diffusion  of  relative  vorticity  with  super-slip  boundary  conditions  applied  at 
the  basin  boundary.  The  second  model,  studied  in  Chapter  3,  is  a  reduced  gravity  quasi- 
geostrophic  model  in  a  basin  with  a  zonal  to  meridional  aspect  ratio  of  9:7.  Like  the  first 
model,  it  has  an  explicit  dissipation  operator  with  bottom  drag  and  biharmonic  diffu¬ 
sion  of  vorticity,  but  the  boundary  condition  is  different.  Free-slip  instead  of  super-slip 
boundary  conditions  are  applied  at  the  basin  walls.  In  addition  to  these  differences,  the 
symmetry  of  the  imposed  mechanical  forcing  is  relaxed,  and  an  additional  parameter  is 
added  to  the  model. 

Both  models  have  previously  been  studied,  but  only  through  prognostic  integration 
of  the  time-dependent  model.  The  major  novelty  of  this  work  with  respect  to  previous 
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results  is  that  pseudo-arclength  continuation,  based  on  Newton’s  method,  is  used  to 
find  both  stable  and  unstable  steady  state  solutions.  Marshall  (1984),  computed  time- 
dependent  solutions  of  the  barotropic  model  in  order  to  study  the  role  played  by  eddies 
in  setting  the  mean  state  of  the  circulation  by  redistributing  vorticity  across  the  inter¬ 
gyre  boundary.  McCalpin  and  Haidvogel  (1996),  computed  time-dependent  solutions  of 
the  reduced  gravity  model  in  order  to  study  the  intrinsic  low-frequency  variability  of  the 
circulation  driven  by  a  steady  wind-stress  curl. 

4.1.1  Results  Concerning  the  Low-frequency  Variability  and 
the  Time-mean 

The  thesis  shows  that  much  of  the  low-frequency  variability  of  the  model  is  associated 
with  trajectories  that  tend  in  the  direction  of  unstable  fixed  points.  In  Chapter  3,  we 
show  that  the  multiple  regimes  identified  in  the  time-dependent  simulation  (McCalpin 
and  Haidvogel  (1996))  have  flow  fields  which  are  very  similar  to  the  steady  state  solutions. 
Averaging  the  flow  field  within  each  regime  shows  that  both  the  meandering  structure 
and  the  eastward  extent  of  the  jet  are  remarkably  similar  to  some  of  the  steady  state 
flow  fields.  Many  aspects  of  the  low-frequency  variability  can  be  captured  by  spatial 
structures  in  phase-space  which  point  away  from  the  time-mean  state  and  towards  the 
model  s  fixed  points  (Chapter  3,  Section  3.4.1).  Four  spatial  structures,  constructed 
to  span  the  hyper-plane  formed  by  taking  the  difference  from  the  time-mean  flow-field 
and  four  of  the  fixed  points  can  capture  30%  ^  of  the  total  interface  height  variance  —  a 
significant  amount  considering  that  the  full  system  is  described  by  24881  modes  or  grid 
points. 

There  is  of  course  no  general  reason  for  an  arbitrary  dynamical  system  to  behave  in 
such  a  manner.  It  is  a  property  of  the  ocean  model  in  the  explored  parameter  range 
that  some  of  the  fixed  points  we  found  are  related  to  the  model’s  attractor.  There  is 

^The  partition  of  the  variance  among  the  four  modes  depends  on  the  order  in  which  they  were 
orthonormalized,  but  the  total  is  independent  of  this  ordering 
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no  guarantee  that  the  fixed  points  will  remain  close  to  the  attractor  in  other  parameter 
ranges.  However,  the  investigated  parameter  range  is  important  since  it  produces  realistic 
oceanic  flow  fields. 

The  bifurcation  structure  we  have  mapped  out  is  also  useful  in  organizing  some  of 
the  time-dependent  behavior  of  the  models  as  the  external  parameters  are  varied.  The 
appearance  and  disappearance  of  steady  state  solutions  via  saddle  node  bifurcations  in 
parameter  space  can  mark  the  transition  between  different  time-dependent  dynamical 
regimes. 

In  Chapter  2,  we  found  that  a  saddle  node  bifurcation  which  marks  the  disappearance 
of  non-symmetric  vorticity  exchanging  steady  solutions  also  marks  a  change  in  the  na¬ 
ture  of  the  time-dependent  solutions.  On  one  side  of  the  saddle  node  bifurcation  where 
the  non-symmetric  steady  solutions  exist,  turbulent  time-dependent  solutions  with  in¬ 
stantaneous  fields  similar  to  the  steady-state  flow  fields,  transfer  vorticity  across  the  line 
of  zero  wind-stress  curl.  This  transfer  allows  a  global  vorticity  balance  which  retains  a 
time-mean  Sverdrup  solution  in  parts  of  the  basin.  On  the  other  side  of  the  saddle  node 
bifurcation  where  the  non-symmetric  solutions  do  not  exist,  all  trajectories  converge  to 
an  anti-symmetric  steady  state  solution  with  basin  filling  inertial  gyres.  In  this  case  the 
Sverdrup  solution  is  completely  destroyed.  This  suggests  that  the  part  of  the  attractor 
that  contains  the  turbulent  vorticity  exchanging  solutions  is  linked  to  the  presence  of  the 
non-symmetric  fixed  point. 

In  Chapter  3,  a  saddle  node  bifurcation  marking  the  disappearance  of  a  high  energy 
nearly  anti-symmetric  steady  state  solution  separates  two  regimes  with  different  low- 
frequency  variability  characteristics.  On  one  side  of  the  saddle-node  bifurcation  point 
where  the  wind-stress  profile  is  more  nearly  anti-symmetric,  the  model  exhibits  low- 
frequency  variability  associated  with  irregular  transitions  to  a  quasi-stable  high  energy 
state.  On  the  other  side  of  the  saddle-node  bifurcation  point,  where  the  wind-stress 
profile  is  more  asymmetric,  no  such  low-frequency  variability  exists.  A  further  result 
is  that  the  closer  the  model  trajectory  approaches  the  fixed  point  the  longer  the  quasi- 
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stable  high  energy  regime  persists.  It  seems  therefore  reasonable  that  the  fixed  point 
attracts  the  time-dependent  trajectories  along  its  stable  manifold  to  produce  the  irregular 
high  energy  episodes  and  the  associated  low-frequency  variability.  As  the  saddle-node 
bifurcation  point  is  approached,  the  stable  manifold  becomes  progressively  less  attracting, 
which  makes  the  high  energy  events  increasingly  less  likely.  Past  the  bifurcation  point, 
the  elimination  of  the  fixed  point  also  eliminates  the  mechanism  that  generates  the  low- 
frequency  variability. 

4.1.2  Results  Concerning  the  Aspect  Ratio  of  the  Basin 

A  typical  property  of  the  bifurcation  structure  for  fluid  dynamics  problems  is  that  it 
depends  strongly  on  the  underlying  geometry  of  the  model  (Seydel  (1994)).  This  is  also 
the  case  in  this  study.  The  different  zonal  to  meridional  aspect  ratio  of  the  basins  used 
in  the  two  models  has  dramatic  effects  on  the  bifurcation  structure  of  the  steady  state 
solutions.  As  the  flow  is  made  more  nonlinear  by  either  decreasing  the  dissipation  or 
increasing  the  forcing,  the  recirculation  cells  expand  in  size  in  the  zonal  direction  until 
they  reach  the  eastern  wall.  Successive  symmetry  breaking  pitchfork  bifurcations  occur 
as  the  jet  progresses  deeper  into  the  basin  interior.  The  larger  zonal  to  meridional  aspect 
ratio  used  in  Chapter  3,  (9:7),  allows  8  symmetry  breaking  pitchfork  bifurcations  to  occur 
before  the  jet  reaches  the  eastern  wall.  The  previous  studies  of  Jiang  et  al  (1995),  Cessi 
and  lerley  (1995),  and  Dijkstra  and  Katsman  (1996)  all  used  basins  with  a  1:2  aspect 
ratio.  Cessi  and  lerley  (1995)  found  only  2  pitchfork  bifurcations  in  a  narrow  basin. 
With  a  similar  aspect  ratio  but  with  different  boundary  conditions  the  model  of  Chapter 
2,  also  had  only  2  or  4  pitchfork  bifurcations  depending  on  the  strength  of  the  bottom 
friction. 

The  wider  basin  not  only  allows  additional  pitchfork  bifurcations  to  occur,  but  also 
clarifies  their  nature.  Each  successive  symmetry  breaking  pitchfork  bifurcation  occurs 
when  the  viscosity  is  reduced  enough  for  the  jet  to  penetrate  sufficiently  far  into  the  basin 
interior  to  allow  an  integral  number  of  half  meanders  in  the  jet  in  the  new  equilibria. 


140 


The  length  of  the  meanders  decreases  downstream.  For  the  model  configuration  used  in 
Chapter  3,  each  successive  bifurcation  occurs  when  the  jet  has  penetrated  approximately 
450  additional  kilometers  into  the  basin  interior. 

The  wider  basin  also  allows  for  a  transition  region  with  damped  Rossby  waves  in 
between  the  recirculation  cells  and  the  Sverdrup  interior.  This  Rossby  wave  field  is 
similar  to  the  solution  proposed  by  Moore  (1963)  for  the  structure  of  the  inertial  western 
boundary  layer  where  the  Sverdrup  flow  is  eastward.  As  discussed  by  Pedlosky  (1996), 
the  Moore  solution  cannot  be  regarded  as  a  model  of  the  western  boundary  layer  but 
should  be  viewed  as  a  distinct  dynamical  regime  for  the  region  separating  the  recirculation 
cells  and  the  Sverdrup  interior.  Before  the  present  study,  the  stationary  wave  field  had 
only  been  observed  in  a  model  with  no-slip  boundary  conditions,  because  only  with  such 
boundary  conditions  does  the  recirculation  cell  remain  limited  enough  to  allow  for  a  zonal 
Sverdrup  interior  that  can  support  stationary  Rossby  waves.  With  a  wider  basin,  the 
stationary  Rossby  wave  field  can  exist  even  with  free  slip  boundary  conditions. 

Finally,  it  is  interesting  to  note  that  the  saddle-node  bifurcation  leading  to  multiple 
anti-symmetric  equilibria  occurs  at  the  same  parameter  value  for  which  the  jet  reaches 
the  eastern  wall.  This  coincidence  is  not  so  striking  for  the  narrow  basin  experiments. 
This  result  suggests  that  the  cusp  catastrophe  identified  by  lerley  and  Sheremet  (1995) 
is  related  to  the  interaction  of  the  jet  with  the  eastern  wall. 

4.2  Future  Work 

The  results  presented  in  this  thesis  have  pointed  out  the  major  role  played  by  the  ge¬ 
ometry  of  the  basin  and  by  the  geometry  of  the  wind-stress  in  modifying  the  bifurcation 
structure  of  the  model.  The  effects  of  non-zonal  winds  and  of  basins  with  different  aspect 
ratios  could  be  easily  analyzed  using  the  present  method.  The  effect  of  bottom  topogra¬ 
phy  could  also  be  studied  with  the  present  method,  provided  the  topographic  slopes  are 
gentle  enough  to  allow  the  quasi-geostrophic  approximation  to  remain  valid.  Further,  the 
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effects  of  irregular  coast-lines  could  be  studied  by  using  a  finite-elements  discretization 
of  the  governing  equations  so  that  the  basin  geometry  could  be  modified  in  a  continuous 
way,  thereby  making  continuation  methods  applicable  to  the  problem.  Other  effects,  like 
stratification,  could  in  principle  also  be  included  provided  enough  computer  memory  is 
available  to  apply  Newton’s  method. 

The  construction  of  a  damped  stationary  Rossby  wave  model  to  explain  the  nature 
of  the  meandering  structure  of  the  null  modes  at  the  pitchfork  bifurcations  would  useful. 
In  addition,  it  would  be  interesting  to  compare  those  modes  with  those  obtained  by 
assuming  that  the  jet  is  a  zonal  parallel  flow. 
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Appendix  A 


Bifurcation  theory  terminology 


If  for  some  parameter  value  one  of  the  eigenvalues  in  the  spectrum  vanishes  (see  Ap¬ 
pendix  B  for  a  definition  of  the  spectrum),  the  steady  state  solution  and  the  corresponging 
value  for  the  parameter  form  a  bifurcation  point,  or  more  precisely  a  zero-eigenvalue 
bifurcation  point.  The  most  common  are  the  saddle-node,  transcritical  and  pitch- 
fork  bifurcations.  Figure  A.l  gives  typical  graphs  for  each  of  these  bifurcation  types 
and  introduces  some  alternative  terminology  that  is  sometimes  used. 

In  a  2-dimensional  parameter  space  there  is  an  unfolding  of  the  fold  bifurcation  at  a 
cusp  point.  This  folding  or  unfolding  is  also  known  as  a  fold  catastrophe.  Figure  A. 2 
shows  how  a  the  region  where  the  surface  folds  over  on  itself  projects  as  a  cusp  on  a 
2-demensional  parameter  space. 

The  symmetry  breaking  pitchfork  bifurcation  can  also  be  understood  in  terms  of 
similar  folded  surface  in  3-dimensions.  Pitchfork  bifurcations  occur  when  there  is  a 
symmetry  in  the  problem.  Suppose  that  a  is  the  parameter  controling  the  anti-symmetry 
of  the  problem  such  that  for  0  =  0,  the  problem  is  perfectly  anti-symmetric,  and  for  a  0 
the  symmetry  is  destroyed.  The  pitchfork  bifurcation  occurs  on  the  intersection  of  the 
folded  surface  with  the  plane  defined  by  o:  =  0.  The  perturbed  or  imperfect  pitchfork 
bifurcation  appears  on  the  intersection  of  the  folded  surface  with  a  plane  define  by 
a  =  const  7^  0.  Figure  A. 3  shows  the  pitchfork  bifurcation  when  a  =  0  (perfect 
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Figure  A.l:  Typical  zero-eigenvalue  bifurcations.  The  control  paramter  varies  along  the 
abscissa  and  the  solution  varies  along  the  ordinate. 
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Figure  A.3:  Left  panel:  pitchfork  bifurcation  for  a  =  0  (perfect  anti-symmetry).  Right 
panel:  perturbed  pitchfork  bifurcation  for  a  >  0  (no  symmetry). 

anti-symmetry)  and  the  perturbed  pitchfork  bifurcation  for  a  >  0  (no  symmetry). 

Finaly,  if  a  pair  of  complex  conjugate  eigenvalues  crosses  the  imaginary  axis,  and  all 
other  eigenvalues  are  stable  (i.e.  to  the  left  of  the  imaginary  axis),  a  Hopf  bifurcation 
occurs.  The  Hopf  bifurcation  leads  to  self-sustained  oscillations. 


Appendix  B 


Method  of  solution 


The  governing  equation  is  discretized  using  finite  differences  on  a  nonuniform  rectangular 
grid,  with  grid  points  in  the  x  direction  and  Ny  grid  points  in  the  y  direction.  The 
scheme  is  made  to  be  second  order  accurate  on  the  nonuniform  grid  by  using  a  method 
outlined  in  Marti  et  al.  (Marti  et  al.  1992).  The  physical  coordinates,  {x,y)  are  given  in 
terms  of  the  computational  coordinates  {(*,  j)  \  i  =  1,2, ...  Ny  and  j  =  1,2,...  N^}  by 
the  following  formulae 


X{jf-m,X{j)  ^  Y{i)^  +  myY(i) 

X  = - 7— - ,  and  y  =  -  - 


1  + 


1  + 


where 


XU) 


N,-l 


N,-l 


(B.l) 


(B.2) 


The  quantities  nix  and  ruy  are  adjustable  parameters  that  control  the  degree  to  which  grid 
lines  are  concentrated  near  the  western  boundary  and  near  the  line  of  zero  wind-stress 
curl.  See  Figure  B.l  for  an  example  of  a  coarse  computational  grid  {tUx  =  niy  =  1.5) 
used  in  Chapter  3. 

The  Jacobian  is  discretized  using  Arakawa’s  (Arakawa  1966),  formulation  with  the 
appropriate  modification  for  the  nonuniform  grid  spacing  (Salmon  and  Talley  1988). 
For  the  super-slip  boundary  condition  used  in  Chapter  2,  the  potential  vorticity  on  the 
boundary  is  treated  as  an  unknown  that  must  be  solved  as  part  of  the  solution,  and  the 
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Coarse  grid 


Figure  B.l;  Computational  grid  used  for  coarse  resolution  solutions.  The  grid  spacing 
in  the  X  direction  varies  from  30  km  along  the  western  wall  to  90  km  along  the  eastern 
wall,  and  the  grid  spaceing  in  the  Y  direction  varies  from  30  km  in  the  center  of  the  basin 
to  90  km  along  the  northern  and  southern  walls. 

no  flux  boundary  conditions  are  treated  using  second  order  accurate  centered  differences. 
After  discretization,  the  PDF  is  expressed  as  a  coupled  system  of  nonlinear  ordinary 
differential  equations  for  the  time  dependent  case  and  a  coupled  system  of  nonlinear 
algebraic  equations  for  the  steady  case.  There  is  one  unknown  for  each  grid  point  which 
can  be  organized  into  a  state  vector  u.  Elements  of  u  corresponding  to  grid  points  in  the 
interior  of  the  domain  are  the  values  of  il)  evaluated  at  the  interior  grid  points  and  those 
elements  of  u  corresponding  to  grid-points  on  the  boundary  of  the  domain  are  the  values 
of  C  evaluated  at  the  boundary  grid  points.  The  discretized  equation  for  the  steady  state 
can  be  written  as  follows 

F{u-,  Si,  6s,  6h)  =  0.  (B.3) 

The  steady  state  solutions  are  found  using  a  pseudo-arclength  predictor-corrector 
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continuation  algorithm.  An  Euler  predictor  step  is  used  as  an  initial  estimate  of  the  solu¬ 
tion,  followed  by  an  iterative  Newton  corrector.  The  solution  branches  are  parameterized 
using  pseudo-arclength  continuation  (Seydel  1994),  which  allows  the  solution  branches 
to  be  traced  past  singular  points.  The  pseudo-arclength,  s,  is  given  for  0  <  C  <  1,  by  the 
folowing  parameterizating  equation 

N 

0  =  p(u;  5h,  s)  :=  (s  -  sj)  -  (1  -  -  6H{sj))dui{sj)/ds  -  C  -  Ui{sj))dui(sj)/ds. 

(B.4) 

In  the  above  equation,  Sh  can  be  substituted  by  either  5i,  5s,  or  any  other  parameter. 

Biffurcation  points  leading  to  multiple  equilibria  occur  at  parameter  values  where  one 
of  the  eigenmodes  of  the  system  linearized  about  the  fixed  point  is  zero.  In  other  words, 
if  we  substitute  solutions  of  the  form 

i’'ix,y,t)  =  e^^<f>{x,y),  (B.5) 

into  the  linearized  system  of  equations  we  obtain  the  following  eigenvalue  problem 

aVV  +  JO’S,  5]V^<I>)  +  J{^,  +  y)  =  -SsV^<f>  -  5^V"(VV),  (B.6) 

where  tps  is  the  steady  state  equilibria.  The  set  of  eigenvalues  {a}  is  called  the  spectrum. 

The  discretizated  linear  stability  equation  can  be  rewritten  in  matrix  form  as  a  gen¬ 
eralized  eigenvalue  problem 

FuV  =  crLv.  (B.7) 

Fu  is  the  Jacobian  matrix  of  partial  derivatives  and  L  is  a  discretized  version  of  the 
Laplacian  operator.  The  discretized  eigenmode  is  given  by  v,  and  its  eigenvalue  is  given 
by  a. 

Zero-eigenvalue  bifurcations  can  be  detected  without  solving  the  stability  problem. 
Instead,  the  sign  of  the  determinant  of  the  Jacobian  matrix  of  partial  derivatives  of  the 
discretized  system  of  equations  is  monitered  -  any  change  of  sign  indicates  that  a  real 
eigen  value  has  crossed  the  imaginary  axis. 


149 


To  detect  if  a  Hopf  bifurcation  as  occurred  as  one  of  the  equation’ s  parameters  is 
varied,  a  method  introduced  by  Neubert  (Neubert  1993),  is  used.  The  method  consists 
of  a  predictor-corrector  strategy  to  follow  the  curve  of  the  dominant  eigenvalue  o  as 
a  function  of  one  of  the  parameters.  A  subsequent  computation  is  used  to  detect  the 
possibility  of  the  occurrence  of  an  exchange  of  roles  whereby  an  eigenvalue  not  being 
followed  becomes  dominant  with  respect  to  the  real  part. 

B.l  Some  practical  details 

Implementing  the  above  procedure  -  i.e.  setting  up  the  operator  Fu  -  can  be  greatly 
simplified  by  (i)  use  of  a  modular  coding  approach  in  setting  up  the  operators  and  (ii) 
use  of  MATLAB. 

Modular  coding  of  linearized  operators 

The  Fu  matrix  is  rather  complex  —  but  the  following  modular  approach  brings  some 
elegance  to  bear  on  the  problem  and  should  greatly  reduces  the  chance  of  coding  errors. 

Let  us  set  ou  the  state  vector  ordered  in  the  following  way,  in  a  great  big  column 
vector  of  dimension  {I  x  N^Ny) 
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Define:  the  ‘forward  differencing  operator’,  an  (n  x  n)  matrix: 


the  ‘backward  differencing  operator’ 


the  ‘Dirichelet  b.c.  operator’ 


the  ‘periodic  difference  operator’ 


Then  the  centered-difference  second  derivative  operator  with  Dirichelet  b.c.  can  be  writ¬ 
ten: 

The  centered-difference  second  derivative  operator  with  periodic  b.c: 

=  +■*»]) 

Laplacian  operator  in  a  channel  with  Dirichelet  b.c.  on  the  northern  and  southern  walls 


Note  that  A  is  an  (n^nj,  x  UxTiy)  matrix  and  and  {ux  x  Ux)  matrix  etc. 

An  example 

Centered  difference  advection  in  the  x-direction,  by  a  constant  zonal  flow  U: 

^  ^  ^ny  -  ®  [P+^  +  (5+J  +  [P^"  -h  =  A 


152 


where  U  is  an  {rixTiy  x  UxTiy)  matrix. 


Advection  of  vorticity,  for  example,  can  be  represented  thus: 

OX 

and  likewise  for  6tc  etc. 


Appendix  C 


Convergence  and  grid  resolution 


In  this  section,  the  issue  of  convergence  of  the  steady  state  solutions  as  the  grid  resolution 
is  increased  is  discussed.  Some  of  the  computations  were  repeated  using  a  uniform  grid 
with  half  and  twice  as  many  grid  points  in  both  the  x  and  y  direction.  The  three  grids 
had  (33  x  17),  (65  x  33),  and  (129  x  65)  grid  points  in  the  x  and  y  directions.  In  all  the 
computations  5s  and  5}j  are  fixed  at  0.01  and  0.04  respectively. 

The  bifurcation  structure  found  in  Chapter  2  is  not  changed  qualitatively  in  going 
from  a  33  X  17  grid  point  model  to  a  129  x  65  grid  point  model.  The  coarse  resolution 
model  with  only  33  x  17  grid  points  is  sufficient  to  capture  the  existence  of  all  the 
bifurcations  of  the  anti-symmetric  branch  found  using  the  higher  resolution  model  with 
129  X  65  grid  points.  Table  1  shows  that  the  discretized  model  appears  to  be  converging 
quantitatively  as  the  grid  spacing  is  reduced. 
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Bifurcation  Point 

Resolution 

Si 

0{}P) 

relative  difference 

PF  N1 

33  X  17 

0.0362 

2% 

65  X  33 

0.0334 

0.5% 

129  X  65 

0.0340 

PF  N2 

33  X  17 

0.0743 

0.2% 

65  X  33 

0.0737 

0.1% 

129  X  65 

0.0734 

Hopf  1 

33  X  17 

0.0762 

0.6% 

65  X  33 

0.0743 

0.1% 

129  X  65 

0.0740 

Hopf  2 

33  X  17 

0.1720 

3% 

65  X  33 

0.1522 

1% 

129  X  65 

0.1514 

PF  N3 

33  X  17 

0.1416 

8% 

65  X  33 

0.0991 

<0.1% 

129  X  65 

0.0990 

PF  N3 

33  X  17 

0.1577 

7% 

65  X  33 

0.1145 

0.1% 

129  X  65 

0.1140 

Table  C.l:  Comparison  of  the  location  of  the  bifurcation  points  for  6s  =  0.01  and 
6h  =  0.04,  computed  on  three  grids  with  uniform  grid  point  spacing  and  with  33  x  17, 
65  X  33,  and  129  x  65  grid  points  in  the  y  and  x  directions.  Column  3  gives  0{h‘^)  estimate 
for  the  location  of  the  bifurcation  points.  Columns  4  gives  the  relative  difference  of  the 
location  of  the  bifurcation  points  calculated  on  the  different  grids. 
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